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PROGRAMMER'S MANUAL FOR MMLE3, 


A GENERAL FORTRAN PROGRAM FOR 
MAXIMUM LIKELIHOOD PARAMETER ESTIMATION 


Richard E . Maine 
Dryden Flight Research Center 


INTRODUCTION 


This report is a programmer's manual for the FORTRAN IV computer program 
MMLE3, a maximum likelihood parameter estimation program capable of handling 
general bilinear dynamic equations of arbitrary order with measurement noise and/ 
or state noise (process noise) . The basic MMLE3 program is quite general and , 
therefore , applicable to a wide variety of problems . The basic program can interact 
with a set of user-written problem-specific routines to simplify the use of the program 
on specific systems . A set of user routines for the aircraft stability and control 
derivative estimation problem is provided with the program . A companion document , 
the User's Manual (ref. 1) , describes the theory and use of the program. This 
paper contains program listings and suggestions for implementation on various com- 
puter systems . Enough information is given about the purpose and operation of each 
subroutine so that users can make modifications if desired . Complete listings and 
reference maps of the routines are included on microfiche as supplement 1. Four test 
cases are discussed; listings of the input cards and program output for the test 
cases are included on microfiche as supplement 2 . 

It is advised that sections 1 and 2 of this paper be read carefully before attempting 
to implement the MMLE3 program on a computer system. The remainder of the paper , 
particularly appendix A , is a reference for detailed information about the structure 
and coding of the program . 


NOMENCLATURE 


A 

a 


n 


X 


B 

C 

C 


m 

S 

Cy 

D 

E 

FF* 

GG* 

<S<g* 

H 

K 

n 

P 

Q 


state equation matrix 
normal acceleration , g 

longitudinal acceleration , g 

state equation matrix 
observation matrix 
axial force coefficient 

lift coefficient 

rolling moment coefficient 

pitching moment coefficient 

normal force coefficient 

longitudinal force coefficient 

lateral force coefficient 

observation matrix 
observation matrix 

state noise power spectral density matrix 
residual covariance matrix 
measurement noise covariance matrix 
observation matrix 

Kalman filter gain matrix (program variable name KGAIN) 
state noise vector 
Riccati covariance matrix 
pitch rate, deg/ sec 
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A 
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Superscript: 


i 

I 

2 2 

dynamic pressure , N/m (Ibf/ft ) ! 

I 

state equation matrix 

state equation matrix 

time , sec 

control vector 

velocity , m/sec (ft/sec) 

forcing function in state equation 

forcing function in observation equation 

state vector 

corrected state vector 

observation vector 

! 

predicted observation vector j 

angle of attack , deg i 

angle of sideslip , deg 
time interval , sec 

aileron deflection , deg j 

I 

elevator deflection, deg j 

i 

measurement noise vector 

pitch angle , deg 

bank angle , deg 

integral of the transition matrix 

gradient (row vector) 

transpose 


3 


Subscripts: 

i 


a, a , p, 5^, 6g 


0 


general index 


derivative with respect to indicated quantity, per deg 
^ 2 

or per deg 

derivative with respect to rate of change of angle of 
attack, per rad/sec 

bias 


Prefix to matrix names; 

APR a priori weighting 

Suffixes to matrix names: 


I 

L 

M 

N 


inverse 

dimensionalization addition 
dimensionalization ratio 
nondimensional 


V 

variation 

Computer labels: 


ALPHA 

angle of attack , deg 

ALT 

altitude, m (ft) 

AN 

normal acceleration , g 

AX 

longitudinal acceleration , g 

AY 

lateral acceleration , g 

BETA 

angle of sideslip , deg 

DELTA-A 

aileron deflection , deg 

DELTA-E 

elevator deflection , deg 

DELTA-R 

rudder deflection, deg 

MACH 

Mach number 

P 

roll rate , deg/ sec 


4 



1.0 


PHI 

bank angle , deg 

Q 

pitch rate, deg/ sec 

Q-BAR 

2 2 

dynamic pressure, N/m (Ibf/ft ) 

R 

yaw rate, deg/ sec 

THETA 

pitch angle , deg 

V 

velocity, m/sec (ft/sec) 


1.0 FORMAT OF DECKS, LISTINGS, AND MODIFICATIONS 


In this section, the conventions used for reference to specific cards in the MMLE3 
program are defined . The program is maintained at the Dryden Flight Research 
Center as an UPDATE (ref. 2) file under the CDC SCOPE 3.4 and NOS 1.4 operating 
systems . Users familiar with CDC computers will recognize the UPDATE format and 
the FORTRAN extended reference maps (ref. 3). Supplement 1 is a microfiche list- 
ing of the MMLE3 program and reference maps. Appendix A contains detailed 
discussions of individual subroutines and common blocks . 


1 . 1 Decks and Common Decks 


This section describes the format of the card decks of the MMLE3 program. 

The format of these decks is as an UPDATE (ref. 2) source file. This format is 
used for the relative ease of changing matrix dimensions (see sec. 3.2) . 

Users with access to UPDATE will find it most convenient to maintain the pro- 
gram as an UPDATE library; the card decks provided can be used directly as input 
to UPDATE . For users without access to UPDATE , the program COMSUB is provided 
(appendix B) to translate the UPDATE source decks into FORTRAN code. Most of 
the sophisticated features of UPDATE are avoided, so COMSUB is a very simple 
program . Its only function is common deck substitution . 

A common deck is a group of cards that can be copied into several different 
subroutines . Although the name is similar , it has no direct relation with FORTRAN 
COMMON statements; FORTRAN COMMON statements are one convenient application 
of UPDATE common decks . The main advantage of using common decks is in the 
ease of making program modifications . If a common deck is modified , the modifi- 
cations will automatically apply to every copy of that common deck . Since it is not 
rare for a common deck to be copied in 10 or more subroutines , the work of making 
program modifications can be reduced by an order of magnitude . 

The first section of the MMLE3 cards consists of the common decks. Each common 
deck is preceded by an identifying card in the format 

*COMDECK common-deck-name 
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The MMLE3 program and subroutines are placed after the common decks in the 
order shown in appendix A. Each routine (including the main program) is preceded 
by a card in the format 


*DECK deck-name 

This card can be ignored unless UPDATE is used . Wherever a copy of a common 
deck is desired , there is a card in the format 

*CALL common-deck-name 

The program COMSUB substitutes a copy of the appropriate common deck in place 
of this card . 


1 . 2 Reference Maps 


The subroutines in supplement 1 are followed by reference maps . The reference 
maps list the variable names in alphabetical order and give the line number of each 
use of the variable names . The line numbers used in these reference maps are the 
numbers that appear on the left of the listings every fifth line . These line numbers 
are used only in the reference maps; all other references to individual cards in this 
report will use the card numbers (sec .1.3). A complete description of the information 
in the reference maps is found in reference 3 . 


1 . 3 Card Numbers 


Except for the reference maps , all references to individual cards use the card 
numbers shown on the right of the listings in supplement 1 (columns 73 to 90) . The 
card numbers consist of an ident name and a number; both are necessary to specify 
a card . Most of the cards within a subroutine have the subroutine name as an ident . 

Cards with an ident different from the subroutine name are either part of a 
modification or a common deck . If the first character of the ident is $ , the card is 
part of a modification . There are no modifications in the current listing; this de- 
scription is included to allow for possible future changes (see sec . 1.4). 

If the ident of a card is not the same as the subroutine name and does not begin 
with $ , that card is part of a common deck . Common decks are described in section 1 . ] 
Any reference to a card in a common deck applies to every copy of the common deck . 
All of the *CALL cards calling for common decks are listed separately , after the list- 
ings of the common decks , before the main program . 


1 . 4 Modifications 


The MMLE3 program is designed as a working tool that can be modified to fit 
specific applications , rather than as an inviolate whole . In addition , it is possible 
that modifications will be necessary to correct program bugs (naturally , we hope 
not) . Therefore , this section describes the conventions to be used by any modifica- 
tions; these conventions are a subset of CDC UPDATE (ref. 2) . Users with access to 
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CDC computers will find it convenient to use UPDATE to implement modifications . 
For other users , the descriptions and example below should be adequate to define 
the actions required to implement any modifications . 

Each group of modifications is preceded by a card in the format 

*ID correction- set-name 

This card defines the ident (sec . 1 . 3) to be used as part of the card number for any 
cards added by this correction set. For the MMLE3 program, the convention has 
been established that the first character of every correction set name will be $ . 
Correction set names are limited to nine characters in length . 

Insertion of new cards into the program is defined by a card in the format 

*I card-number 

immediately followed by one or more cards to be inserted . The card number is in 
the format 


ident .number 

and defines the card after which the new cards are to be inserted . If the card number 
describes a card in a common deck , the new cards are to be inserted in every copy 
of the common deck . Cards may be inserted after any card in the program , including 
cards inserted by previous modifications and cards which call for inclusion of common 
decks . Cards to be inserted may include cards that call for inclusion of common decks 

Deletion of cards from the program is defined by a card in the format 

*D card-number 


or 


*D card-number- 1. card-number-2 

The first format describes a single card to be deleted; the second format is used to 
delete all cards from card number 1 to card number 2 inclusive . The *D card may be 
followed by one or more cards to be inserted in place of the deleted card (s) . The 
number of cards inserted in this manner does not have to equal the number of cards 
deleted; it can be larger or smaller . 

Cards beginning with */ are comments and can be ignored. 

The following simple example should help clarify some of the ideas of this 
section . The correction set reads the time history from a tape with the time in total 
milliseconds instead of hours, minutes, seconds, and milliseconds. The dimensions 
are increased to allow tape records up to 150 words long, instead of 100. 
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*ID $LONGTAPE 

*/ READS LONG TAPE RECORDS 
*/ TIME IN TOTAL MS 
*D RECRD . 2 

COMMON /RECRD/ EOFTH ,T (4) .RECORD (150) 

*D READTH.20 

READ (UDATA) ITMS , (RECORD (I) ,1=1 ,NREC) 

CALL IHMSMS (ITMS ,T) 

Note, in particular, that the card RECRD .2 is changed in both copies of the 
common deck , RECRD (in subroutines READTH and THDATA) . The original sub- 
routines READTH and THDATA are in supplement 1 . The routines resulting from the 
above modification are shown below; the listing of THDATA is truncated since there 
are no changes in the latter part of the subroutine . 




SUBROUTINE RF40TH (INSTAT ) 

READTH 

2 

c 



READTH 

3 

c 


READS ONE POINT OF INFLT TIME HISTORY 

READTH 

4 

c 


STANDARD VERSION FQk CARD Ok TAPE INPUT 

READTH 

5 

c 


IN5TAT GIVES INPUT STATUS 

READTH 

t 

c 


D INDICATES FIRST CALL TO READTH FOR THIS CASE 

READTH 

7 

c 


1 INDICATES SEARCHING FOP A START TIME# BUT NOT FIRST CALL 

READTH 

6 

c 


? INDICATES READING DATA 

READTH 

9 

c 


THIS ROUTINE SHOULD NOT ALTER INSTAT 

READTH 

10 

c 



READTH 

11 



COMMON /FILES/ LCARD^ IP0NCH» UPRINT# UDATA»UT1 » UT2# UTHOUTf UWT» UPLOT 

FILES 

2 



INTEGER UCARD»UPUNCH,UPRINT#U0ATA»UTl>UT2f UTHCUT»UNT>UPLOT 

FILES 

3 



COMMON /INOPT/ CARO, TAPE 

INOPT 

2 



LOGICAL CARD, TAPE 

INCPT 

3 



COMMON /TNOkD/ NREC#ZCHAN( 08),UCHAN( UA),EYCHAN( 20) 

INORD 

2 



INTEGER ZCHAN,UCHAN, EXCHAN 

INORD 

3 



COMMON /RECRD/ ECFTH, T ( A ), R ECORD ( 150 ) 

SLONGTAPE 

1 



LOGICAL FGFTH 

RECRD 

3 



INTEGER T 

RECRD 

4 



COMMON /TAPPOS/ lTM,REk 

TAPPOS 

2 



LOGICAL REW 

TAPPOS 

3 

c 



READTH 

17 



IE(CaRO) GO TO 100 

READTH 

18 



IF(REN) REMIND UDATA 

READTH 

19 



RE AO (UDATA) ITM S H R EwCRD ( I) * !•!, NREC ) 

SLOHGTAPE 

2 



CALL 1HMSMS(ITMS,T) 

SLGNGTAPE 

3 



GO TG 500 

READTH 

21 


100 

READ(UCARO,1000) T, ( R ECORO ( I > , I-X, NR EC ) 

READTH 

22 

c 



READTH 

23 


lOCO 

FDRMAT(3I2,I3,1X,7F10,A/(8F10,A) ) 

READTH 

24 


5CC 

RETURN 

READTH 

25 



FND 

READTH 

26 
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vUBhTUTINF 7HDATA 

THOATA 

2 

c 


THDATA 

3 

c 

kcAD? INPUT lIHE Hl<:TCklES, FTNOj AVERAGES OF ALL SIGNALS. 

TMDATA 

4 

c 


THOATA 

K 

c 

tCOPS W1T« T«l,LG«n DEPEND ON RELATIVE HRCER OF Z# AND EXTRA 

TFDATA 

6 

c 

VARUeLtS. SCALES# BIaSFS# CHANNEL NLHbERS# OP AVERA^^ES 

THDATA 

7 

c 

IN THE CuNpnN BLOCKS. 

THOATA 

8 

c 

references TO Z VARIARLES THEN EXTEND TQ U AND EXTRA. 

THOATA 

9 

c 


THOATA 

10 


C IMMOK /AVGtnii/ ZAVG( Ob)#UAVG( 0A)#FXAVC-( 20)#ZSIG( 08)» 

AVGCCM 

2 


USIG( CA),FXSTG( 20)#ZNINH( 0B)>UNINH( 0A)>EXPINN< 20)# 

AVGCCM 

3 


ZMAXM( 08)#UMAXM( C4)#EXMAXM( 20 

AVGCOH 

4 


COMMON /BiLiN/ US EA VG# T IMV AR# Z ( U8)#l( 0A)#EXTRA( 2C)>0NES( 04) 

BlLIN 

2 


logical USFAV6#1IMVAR 

BILIN 

3 


CuMMoN /CuM/ NCASE#NPU#NPTS(lS)#ITMSTSa5) 

CCM 

2 


COMMON /FILES/ IC AkO# 1 PUNCH# UPRI NT# U DATA# UTl# UT2# UTHOUT# UWT » UMLDT 

FILES 

2 


INTEGER U<‘APD#UPUNCH#UFRINT#U0ATA#UT1#UT2#UTHDUT#UWT#UPLCT 

FILES 

3 


COMMON /HPAONG/ TiTLP(20)»A0ATF»ATIME# 

HEADNG 

2 


$IGLAB(2# 08)#XLAB(2# 07) # C TNL A B ( 2 # C4)#EXLAB(2» ?0) 

HEACNG 

3 


COMMON /^NORO/ NREC#ZCHAn( 06)#UCHAN( 04)#EXCHAN( 20) 

INCRO 

2 


INTcGEK ZCHAN,UCHAN#EXCHAN 

INCRD 

3 


COMMON /1NTF6R/ DT»NEAT 

INTF6R 

2 


COMMON /naXI^S/ MAXV#MAXZ#MAXU#MAXd#lFX#LCRD 

MAXIMS 

2 


CnhMHN /MOOCuM/ UMOD 

MODCOM 

2 


LOGICAL UMOr 

MCOCOM 

3 


common /R^CPO/ EOFTH#T(4),KECORO<150) 

SLCNGTAPE 

1 


LOGICAL E^'FTH 

RECRO 

3 


INTEGER T 

RECRD 

4 


CCMMON /TAPFOS/ TTM#RER 

TAPFOS 

2 


LOGICAL RPW 

TAPPOS 

3 


COMMON /THOATA/ STC ( i 5 ) # ETC ( 15 ) # THI N# PR INU # MAXR EC , 

TCCATA 

2 


2BiAS( 0P)#LRIA^( 04)#£XBIAS( 2C)#ZSCALE( 08)#USCALE( 04)# 

TCOATA 

3 


FXSCALC 20) 

TCOATA 

4 


INTEGER THlN,STf»ETC 

TCOATA 

5 


LOGICAL PRINT! 

tcdata 

6 

r 


THDATA 

23 


LOGICAL FIRST 

THOATA 

24 

c 


THOATA 

25 


1 . 5 Tapes 

The following format will be used for tape transmittal of the MMLE3 program 
unless explicitly requested otherwise. The availability of tape copies of the MMLE3 
program can be ascertained by writing the author . 

Tapes are nine-track 800 BPl labeled tapes. The VSN of the tape is MMLE3T . 
The label is American National Standard Institute (ANSI) standard with the name 
MMLE3. All data are ASCII-coded card images. (EBCDIC code is available on re- 
quest.) Each card image is a fixed-length 80-character record. Records are 
blocked in fixed-length blocks of length 1200 characters . Each block contdins 
exactly 15 records with no padding; records do not span blocks. 

Each file of data on a tape is terminated by a card with "END-OF-FILE-nn" in 
the first 14 columns . The nn is replaced by the file number . The remaining 66 col- 
umns of the card are filled with dashes . Actual system end-of-files are not used 
because of possible incompatibility between computers . 

There are 12 files of data on the tape. The first file is the UPDATE source cards 
for the MMLE3 program as described in section 1.1. The second file is the UPDATE 
source cards for the EISPACK routines used by MMLE3. The third file is the CDC 
segmentation directives . The fourth and fifth files are the COMPUN program and 
the template used by COMPUN as described in appendix C . The sixth file is an 
alternate template for COMPUN , which results in a correction set for CDC UPDATE 
instead of a complete set of common decks . The seventh file is the COMSUB program 
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described in appendix B . The eighth and ninth files contain the program modification 
and input cards for the one -dimensional check case listed in appendix D . The 10th 
and 11th files contain the program modification and input cards for the two longitu- 
dinal aircraft check cases described in appendix D . The 12th file contains the input 
cards for the lateral-directional standard aircraft test case described in appendix D . 

A microfiche listing of the tape contents will be supplied with the tape . The 
total length of the tape is approximately 10,400 card images. 

The following short program can be used to pick a file from the tape . The file 
number desired is specified on a single card in 110 format. The tape is assigned 
FORTRAN unit number 11 and the file requested is written to FORTRAN unit num- 
ber 12. 


PROGRAM PICK< INPUT* OOTPU1,MMLE3T»FILEN» 

TAPei*iNPUl,TAPe3«CCTPU1*TAPEi.u-MMLf 3I,T*PEli-FUEN) 

C 

C PICK A "FIlc" off of MMLE3I TAPE AND COPY IT TO FILFN. 

C INPUT IS A SINGLE CARL KITH THE OESIREO flLF NUHRER IN FCRMAT 110. 

r 

INTfGtR UREA0fUPRiNT»llN,UOUT 
REAL CARD(£0)#FN0(3) 

DATA ENu/4HEN0->4HQF-F#4HiLE-/ 

Data UREAU/J/#UPPiNT/3/>UiN/i0/» UOUT/11/ 

C 

c ^^t********^****’^**^ PEaO requested File number# 

READ (UREA0>610W) NRICK 
WRITE (UPRlNT#e300) NPICK 
C 0^^t***^m**^^*0^^^** SKIP PRtCEbDlNG PILES* 

REWIND tUN 
NSKIP • NPICK-l 
IF (NSKIP,LE*0) GO TO 3C0 
jU Zfyj IFILE ■ I » NSKIR 

DO ICARO • I * 20000 

READ (UIN#8900) ( C AKD ( 1 ) > T« 1* 4 ) 

IF (CARD(I)*E0«t:h0U) *AN0* C ARD ( 2 ) • E Q* tND( 2 ) *AND* 
CAR0(3)*EQ«tND(3)) GO TO 
100 CONTINUE 

150 WRITF (URRINT^bJwi) CAK0(4)> iCARD 

201 CONTINUE 

C COPY LtSIREO FILE. 

an 400 icard ■ i > 2 ocuo 

READ (UlN#b90u) CARD 

IF (CARDU )*EQ*END(1) *AND. C ARD ( 2 ) * bO* lNO ( 2 ) *ANn« 

- CARD(3)*EQ*kN0(3)) GO TO 5*^0 

WRITE (U0UT#69f>9) CARD 
400 CONTINUE 

WRITE (UPRINT#83g2) C AR D ( 4 ) , I CAR D 
C 

4100 FORhAT(Ilv) 

83Aa F0RHAT(**0FILE NUN8ER”>13>« REQUESTED PROP THE TAPE.**) 

830t FORHATC FILE •*fA2.*' ^iKI PP cO. **# 16# •* CAROS IN ^UE.**) 

8302 F0RHAT<"0FUE -#A2#** CUPIEO* " # 1 6# *• CAkDS IN PILE.**) 

8VO FuRHAT(2^A4) 

STUP 

END 
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2.0 IMPLEMENTATION CONSIDERATIONS 


This section discusses the considerations in implementing the MMLE3 program 
on various computer systems . The program was designed with generality in mind , 
so code peculiar to specific systems was avoided. The MMLE3 program has been 
checked out on CDC computers with SCOPE 3 .4 and NOS 1.4 operating systems and 
on IBM computers with the OS/360 operating system. It has also been run on Univac 
and Harris equipment . It should not be difficult to implement on any large computer 
system with a FORTRAN IV compiler and a CalComp plotter (ref. 4) . 


2.1 FORTRAN 

The MMLE3 program uses mostly ANSI standard FORTRAN IV (ref. 5) . All 
character data are stored four characters per word for machine generality . The few 
statements used that do not conform to ANSI standards are described below . Compilers 
on most large computer systems accept these usages . 

The program card , which is nonstandard , is necessary on CDC systems to 
define the files to be used . This card appears only in the main program and should 
be deleted if it is not appropriate to a particular system . 

Nonstandard subscripts are used in several places . Most large system com- 
pilers accept the forms used . If a particular compiler does not accept them , the 
nonstandard subscripts can be replaced by dummy variables , defined immediately 
before their use . 

Unsubscripted array names are used in data statements throughout the program. 
This usage, although nonstandard, is very common. The equivalent ANSI standard 
forms are quite tedious but can be substituted if necessary . 

ANSI standard FORTRAN restricts the assignment of values to formal parameters 
in a subroutine if two or more formal parameters are associated with the same actual 
parameter. The MMLE3 program does not abide by this restriction. The restric- 
tion is a quite conservative way of avoiding potential problems in optimizing com- 
pilers. The specific usage in MMLE3 is not prone to such problems. The violation 
of this restriction in the ANSI standard cannot be detected at compilation time and 
thus does not result in compilation errors . 

Quotation marks are used in many subroutines to delimit Hollerith fields in 
FORMAT statements (never in DATA statements) . Most compilers accept quotation 
marks or some other character in this role . A simple character translation program 
can be used to change the quotation marks to acceptable characters if necessary . 

The use of NAMELIST input is the most serious potential problem of FORTRAN 
compatibility for MMLE3. Although NAMELIST input is nonstandard, it is available 
on many computer systems . If NAMELIST is not available , a substitute type of input 
will be necessary . If such a substitute does not embody many of the general features 
of NAMELIST, the program will become very cumbersome to use. In particular, 
because of the large numbers of variables and options available , it is desirable to 
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have default values that are used if no values are input . The user then needs to 
set only the values that are different from the defaults . It is also desirable to have 
the input variables identified by name instead of position , because variable names 
are much easier to remember than positions assigned on the input cards . 

The program does not depend on memory being initialized to 0. However, matrix 
elements which might not be defined are tested to see if they contain the special flag 
value "TEST" (see sec. 3.1) . Therefore, memory must not be initialized to infinite 
or indefinite values on computers that have such values. 


2.2 Files 


Input and output conventions vary from system to system; therefore, some 
flexibility must be built into any program that will be used on different computer 
systems . All input or output in MMLE3 uses the variables in common block FILES 
for unit numbers . The values of these variables are set in subroutine VARDEF . 

The unit numbers can be changed as desired for compatibility with various operating 
systems; for CDC systems, it will be necessary to change the file names in the pro- 
gram card to correspond to the unit numbers used . 

Descriptions of the files used by the program are given in reference 1 , section 3.2. 
On IBM systems , the following DD cards are samples of those necessary to define 
the files . It is assumed that a cataloged procedure is used that defines the card 
reader , card punch , and line printer files . The file numbers and record lengths 
in these samples correspond to the values currently used in the program . 


//GO.FT02F001 DD SYSOUT=B 


//GO.FT04F001 DD DISP=OLD,UNIT=2314,VOL=SER=volume, 

// DCB=(RECFM=VSB,LRECL=420,BLKSIZE=4204) ,DSN=name 

//GO.FT07F001 DD DISP=NEW,UNIT=SYSDA,SPACE=(CYL, (10,2)) , 
// DCB=(RECFM=VSB,LRECL=148,BLKSIZE=4444) ,DSN=INTERNAL 

//GO.FT08F001 DD DISP=NEW,UNIT=SYSDA,SPACE=(CYL, (10,2)) , 
// DCB=(RECFM=VSB ,LRECL=192 ,BLKSIZE=3844) ,DSN=TOPLOT 

//GO.FT09F001 DD DUMMY 
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or 


//GO.FT09F001 DD DISP=(NEW,CATLG) ,UNIT=2314,VOL=SER=volume, 
// SPACE=(CYL,(5,1),RLSE), 

// DGB=,(RECFM=VSB,LRECL=148,BLKSIZE=2964) ,DSN=TH0UT 


/ /GO .FTlOFOOl DD DISP= (NEW ,CATLG) ,UNIT=2314 ,VOL=SER=volume , 
// SPACE=(CYL,(5,1)3LSE). 

// DCB=(RECFM=FB ,LRECL=80 ,BLKSIZE=3200) ,DSN=WTDATA 
(DISP can also be OLD) 


//GO.FT13F001 DD DISP=(NEW,CATLG),UNIT=2314,VOL=SER=volume, 

// SPACE= (CYL ,(5,1) ,RLSE) ,DSN=PLOT 

(on some systems replace FT13F001 by PLOTTAPE) 


On CDC systems , no special control cards are needed except those to attach , request , 
or catalog any permanent files or tapes . 

Core space can be conserved by lowering the I/O buffer sizes . The program 
uses the scratch files UTl and UT2 (7 and 8, respectively) extensively. If the 
buffer size for these two files is lowered too much , execution time will increase . The 
buffer size of file UDATA (4) can also affect execution time if this file is long. 

The remaining files are used seldom enough that even very small buffer sizes will 
not affect the execution time . On CDC computers , buffer size is specified directly 
on the program card . On IBM computers , the total buffer size is determined by 
the block size and number of buffers specified on the DD cards . 


2 . 3 Segmentation or Overlay 

Efficient operation of the MMLE3 program requires the use of segmented or 
overlaid loading. The program requires over 200,000g words of core to execute on a 

CDC computer if segmentation is not used. The use of segmented loading reduces the 
core requirements to approximately 54 , OOOg words . This large difference in core re- 
quirements results because the program was written to take advantage of segmentation . 

The segmentation structure for MMLE3 is shown in figure 1. Each box repre- 
sents a segment . The name assigned to the segment (usually the name of one of the 
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subroutines in the segment) is underlined at the top of the box . A brief description 
of the major functions accomplished in the segment follows. 



Figure 1. Segmentation structure. 
14 
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The segmentation directives required to implement this structure on a CDC 
SCOPE 3 . 4 system are given below . The TREE and LEVEL directives define the basic 
structure , The INCLUDE directives assign subroutines to the appropriate segments , 
and the GLOBAL directives assign the common blocks . All of the subroutines and 
common blocks of the MMLE3 program itself are explicitly assigned by these directives. 
System routines and common blocks are not included, except for the CalComp plotter 
routines and the NAMELIST input routine NAMIN= . The name of the NAMELIST input 
routine will be different on different systems . On some systems , additional directives 
may be required to assign some of the system routines and common blocks to the 
root segment . 


* xriLgl SEGMENTATICN directives 31 JUL 60 RICH MAINE 

• 

TREE MHLE3- (GETSET» NOPLOTS, THPIOT I 
MnPLHTS tree GET-aNTREEfOOTREEl 

INTRPE TREE TITL ES- ( EDIT# HATS ET ) 

OQTRcc tree NEWTON- (GTREE,N6TREE) 

CTREE TREE GIRL-* KALHAN.REAT) 

NGTRFE tree notgirl-(upoate,suhouti 

• 


NNl 

CFTSCT 

GET 

GET 

POIT 

HATSET 

H4TSET 

MATS'^T 

K Af *<AN 

REAT 

RPAT 

NOTGTRl 

DPOATE 

UPOATC 

THPL DT 


INCLUDE MHLE3,AeEN0f MIL 
TNCLt'OE 6ETSET»VARDEFf WTIN 

INCLUDE GET, GETLAB#HE AO, LOCATE, SET #SET1»SET2»SP IT, SPITEP, TOT 1 
INCLUDE Z0T2 

INCLUDE EDIT, DIGIT, USERIN 

TNCl UOF HATSFT..4LL0N, AyEPAG,C0PPAT,CCMN;CGNSTR.GVALVE,HARDC 
INCLUDE INTE RP, LO ADE n,HAT0EF,HAU0,MATNa,MTL0A0, MTS ET, ONCE 
INCLUDE READTH,SETCON,THOATA,THMOD, VARY,WIDEF,WTTRAN 
INCLUDE KALMAN,GRA0K,GRA0P,LYAPCB,R1CATC 

INCLUDE REAT,CALLAM,0IHl,DIH2,EAT,GRAD,GRA0IC,XNIT,aBSEPV 
INCLUDE SPIDIM,VHADO,MAKFL,MAKFH,HAKEVW,THOUT,UINIT 
INCLUDE REDUCE, SINVrSSIMFQ 

INCLUDE UPOATe,APRAOn,BlAS,DFACT,FAOJ,FLlHlT,HOVE,HVHULT 
INCLUDE RESIDS 

INCLUDE $UN0UT,Cff AMER,CRSeT,ERRTHP,OUTPUN,PLOP 
INCLUDE THPLCT, IHHSMS 


GET 

GET 

GET 

GET 

TITLES 

NATSET 

NEWTON 

NEWTON 

GIRL 

GIRL 

KALMAN 

SUHOUT 

♦ FDR 


GLOBAL BUF, COM, 01 MALL, FILES, HE ADNG,INTEGR,MAXCOH,HAXIH, MAXIMS 
GLOBAL MODCQM,OUTOPT,SIZE,TAPPLS,TCPLOT 
GLOBAL FLC0ND,GF0FFS 

GLOBAL AVGCOP,BlLlN,OETERM,DlHhAT, ECGM, FC UM, GICOM, 6R SI ZE 
global lCOND,MAPCOM,MATRAT,HAIRlX,OaSRV,PHlCOH,SQFCtJM, SUMSAV 
GLOBAL TOGIPL 

GLOBAL INERTS, INSTR,LCNLAT 

GLOBAL lNMAT,INOPT,lNnRO,HATIN,MArLAB,TODATA,UVCOM 
GLOBAL RECRO 

GLOBAL ANCOM,OUMCUM,OUNVEC,ERL1ST,GRADS,GROCOH,KCOM,P8COM 

GLOBAL SUMCOM,XSUMS 

GLOBAL TOGRAO,BIASES 

GLOBAL GRAV 

GLOBAL RICCOM 

GLOBAL CRMAT 

GCRATCH STORAGE AS NEEDED 
GLOBAL GLOBAL 


LEVEL 


TREE NAMIN 

NAMTN INCLUDE NANIN-, IOIGIT,ROWCDL 


TREE ADO 

4 on INCLUDE ADD, AOOPAR,OHULT,GETP, GETP AR, IDFNTl, I NV, MAKE, MULT 

Ann INCLUDE MULTT,SHULT, SUB, SUHUL T,SYM, TR4NSP,UNSET, ZHUIT, ZOT 

♦ 

TREE AXES 

Axes INCLUDE AXES,LlNES#PLTOAT,SCALE?,SYMfiLA,TtTPLT 

axes include factor, NUMBER, symbol 

• 

TREE EIG6NG“<ELMHES,H0R,HQR2) 

FLAMES INCLUDE ELMHES, BaLA NC , EL TR AN 

MOR INCLUDE HQR 

NOR? INCLUDE M0R2,BALBAK 

4 


FHD 
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The statement GLOBAL GLOBAL in the above directives declares that the common 
block name GLOBAL will be placed in the root segment so that it can be accessed 
from any subroutine . The program as published does not have a common block 
named GLOBAL . This statement is provided in case the user needs an additional 
common block for some temporary modification . The user can define the common 
block GLOBAL as needed , and the above directives will insure its proper placement 
in the segmentation structure . 

The following directives implement essentially the same structure on the CDC 
NOS 1 . 4 operating system . The routines not mentioned explicitly and all of the 
common blocks are automatically assigned by the system . 


» 11LE3 SEShENTATION OIKECTIVES FOR NOS. 12 SEFT 80 RICH MAINE. 

« 

CTHHQN 

HMLE3-(GETSET»NOPLOTS»THf'LOT) 

H0»>nTS tree GET-(INTR€E#00TREE) 

TNTREE TREE TITLE$-(E0IT#HATSET) 

D0TR*^E TREE NE WTON- ( GTRE E* NGTREE ) 

GTREE TREE GIRL-( KALHAN# RE ATI 

MGTREP TREE SINV-( UPOATE»StlMOUT) 

♦ 

GETSPT include varoee»wtin 

GET INCLUDE GETL AB# SETl# SETZ# SPITt ZOTl^ZOTZ 

M ATGET TNCLUDE QNC E, HTLO AO# TH0ATA> MAT CE F, COPP AT, ALLOW 
UPDATE INCLUDE APR ADO# 0l AS # OFACT# F LIM IT» M VMLLT#R ES ID S 
X AL»*AN INCLUDE 6RA0K 

RE AT INCLUDE 0IM2#GRA0#GR AOICf INIT# 00 SERV # SP ID IH# VMaDO# M AKE VW> THOUT 

♦ 

LEVEL 

* 

TREE NAMIN« 

NAMIN. INCLUDE IDIGITt ROWCOL 

* 

TREE ADO 

4 no INCLUDE aoo#adopar#dmult#getp#6etpar# iuenti#inv# make# mult 

ADD tmcLUDE MULTl#SMULT#SU0#SUMULT#SrM#TRANSP#UNSET#ZMULT#ZOT 

A 

TREE AXES 

Axes INCLUDE AXES #L IMPS# PLTf>AT#SCALE2#SYMBLA#TlTPLT 

AXES INCLUDE FAC TOR# NUMBE R# SYMBOL 

♦ 

TPfcE £IGENG-(ELMHFS#HQR#H0R2) 

EL^Mcs INCLUDE 8 AL ANC # ELTR AN 

NOP2 INCLUDE 8ALBAK 

«N0 


If MMLE3 is to be run without segmentation or overlay, several small changes 
will reduce the amount of core required. First, subroutine WTIN should be removed 
from the MMLE3 program and run as a separate job; alternately, the dimensions in 
WTIN could be decreased . Second , the dimensions in subroutine THPLOT could be 
decreased , at the expense of an increased number of disk reads . Third , some of 
the large common blocks from different segments could be equivalenced . Fourth, 
subroutine HQR2 could be modified to include HQR as a subset. 


2.4 EISPACK Routines 


The MMLE3 program uses the subroutines described in reference 6 for obtaining 
the eigenvalues and eigenvectors of matrices . This subroutine package (EISPACK) 
is widely used and well documented. The EISPACK routines used by MMLE3 are 
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BALANC , BALBAK , ELMHES , ELTRAN , HQR , and HQR2 from release 2 of EISPACK . 
Listings of these EISPACK routines are included in supplement 1, in case the user’s 
installation does not have the EISPACK release 2 library available . The listings 
include brief descriptions of the functions and usage of the routines . 

The EISPACK routine BALANC contains a machine-dependent variable RADIX . 
This variable should be defined at card 64 to be the base of the machine floating 
point representation. For CDC machines, this value is 2 (as shown in the listing 
in supplement 1) . On IBM machines , RADIX should be 16 . RADIX is used to insure 
that the arithmetic in subroutine BALANC will be exact, with no rounding error. 

Subroutines HQR and HQR2 contain a machine-dependent variable MACHEP . 

This variable should be defined at cards HQR. 63 and HQR2.87 to be the smallest 
positive , floating point number which , when added to 1 , gives a result not equal 

-47 

to 1 . For CDC machines this value is 2 , as shown in the listings in supplement 1 . 

Different values would be used for other machines . MACHEP is used in convergence 
tests . 


2 . 5 Plotting 

High-resolution time history plots are needed to adequately evaluate the results 
from MMLE3. Line printer plots are usually inadequate. The MMLE3 program uses 
CalComp plotter software (ref. 4) . If a CalComp plotter is not available, the routines 
calling the CalComp software must be changed. The only CalComp subroutines 
called are PLOTS , PLOT , FACTOR , SYMBOL , and NUMBER . The calls to the CalComp 
software are in the main program , subroutine THPLOT and subroutines AXES , LINES , 
PLTDAT , and SYMBL4 called by THPLOT . 

If the plotting subroutines are rewritten, it is important to recall that MMLE3 
stores labels four characters per word for machine generality . Therefore , each 
word must be treated separately . Subroutine SYMBL4 is an example of this type of 
treatment . On a computer that stores four characters per word , SYMBL4 is not 
needed because a direct call to SYMBOL will work as well . However , the direct 
call to SYMBOL would not be transportable to machines with more than four charac- 
ters per word . 

The functions of the CalComp subroutines will be described here so that the 
user can adapt the program to other plotting software . 


2.5.1 PLOTS (BUF, NBUF, UPLOT) 

Subroutine PLOTS initializes the plotting software . BUF is a scratch storage 
area used by the plotting software , and NBUF is the length of the BUF vector in 
words . Some CalComp implementations do not require a user-allocated buffer , in 
which case BUF need not be dimensioned . UPLOT is the FORTRAN device number 
assigned for the plot file . Some systems ignore UPLOT and assign a system-specified 
file name . 
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2.5.2 PLOT (X, Y, IPEN) 

Subroutine PLOT moves the plotter pen to the position (X , Y) referenced to the 
current origin. IPEN should be ±2, ±3, or 999. If the magnitude of IPEN is 2, the 
pen is down during movement , thus drawing a line . If the magnitude of IPEN is 3 , 
the pen is up during movement . If IPEN is negative , the origin is redefined to be 
the pen position after the move . The value of 999 for IPEN is a special call to close 
the plot file . 


2.5.3 FACTOR (FACT) 

Subroutine FACTOR enlarges or reduces subsequent plots. The size of subse- 
quent plots is FACT times the "normal” size. Calls to FACTOR are not cumulative; 
e.g. , a call with FACT=.25 followed by a call with FACT=.5 makes subsequent plots 
one-half of the original size, rather than one-eighth. The program assumes that the 
original plot size takes all values in inches . 


2.5.4 SYMBOL (X, Y, HOT, I, ANGLE, N) 

Subroutine SYMBOL draws text or symbols . There are two branches in sub- 
routine SYMBOL , depending on the last argument , N . 

If N is positive, N characters taken from the vector I (in A format) are drawn. 

X and Y are the starting coordinates of the lower left-hand corner of the first 
character, and ANGLE is the angle at which they will be drawn. If X or Y is 999, 
the characters are drawn starting at the previous pen location in the corresponding 
axis . HGT is the height of the characters drawn . 

If N is negative , a single symbol is drawn , specified by the integer I . The 
list of symbols is given in reference 4 . The interpretation of X , Y , HGT , and ANGLE 
is the same as when N is positive . If N is -1 , the pen is up during the move to loca- 
tion (X , Y) . If N is -2 or less , the pen is down . The symbols for I values from 0 
to 13 are drawn centered at the location (X, Y); all other symbols treat (X, Y) as 
the lower left-hand corner of the symbol . 

The value N=0 is not used by the MMLE3 program. 


2.5.5 NUMBER (X, Y, HGT, A, ANGLE, N) 

Subroutine NUMBER draws the value of a floating point number . X , Y , HGT , 
and ANGLE are identical to the corresponding arguments in subroutine SYMBOL . A 
is the floating point number to be drawn . N controls the format used . If N is greater 
than or equal to 0 , N digits after the decimal point will be drawn . If N is - 1 , the 
integer part of A is drawn with no decimal point. Values of N less than -1 are not 
used by MMLE3. All numbers drawn are rounded rather than truncated values. 
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2.6 Date and Time 


Calls to the DATE and TIME subroutines are provided to help in identifying 
printed output and plots. Naturally, these subroutines are machine specific. If the 
system does not supply these functions , dummy routines which return blanks can be 
used. The calls to these routines occur in subroutines TITLES and PLTDAT. Any 
other useful identifying information may be substituted in these places , if desired . 


2 . 7 Assembly Language Routines 

The MMLE3 program is coded entirely in FORTRAN . No assembly language 
routines are included in the listings of supplement 1 . The program spends a 
considerable portion of its CPU time in three matrix multiplication routines— MULT , 
SUMULT , and ZMULT . These routines are quite short and simple . It is , therefore , 
worth considering the use of assembly language replacements for these three 
routines in order to decrease the program run time. On a CDC Cyber 70/73, assem- 
bly language (COMPASS) versions of these subroutines run in about one-half of 
the CPU time required for the FORTRAN routines . This results in a decrease in the 
overall program run time of about 20 percent. Assembly language replacements 
for other subroutines are not generally worth the effort . 

Listings of the COMPASS routines are shown below . The macros ENTRYP and 
CALL used in these routines are also listed . 
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jncMT MULT 


CJMKCNT maipix Ml^TIPLlC/iTICN C*A+P 



tSF 

DATA 

MAX 

P- ' 

10 

Al 

0SS 

IP 

JJ 

ess 

10 

MT X 

BS^ 

IP 

J 

PiS 

IP 

KM 

OSS 

10 


ENTRY® MULT 

SAVE CuRHAL 

pakametfrs 


SA2 

A141P 


SA3 

A2 + 1B 


BX^> 

XI 


BX7 

X2 


SA6 

ARGLl 


SA7 

ARGL2 


BX6 

X3 


SA6 

ARGL3 


CAtl 

GET»(Of^AX>lI>JJI, ARGLl 


CALL 

GE I f (0* MIX, j, KK) , ARGL2 


CALL 

SFTif (OfMAXf ITfKK), ARGL3 

LGHf AF t 

VALUFS OF 4 ANO JJ 


SA? 

J 


5AA 

JJ 


IXO 

X<»-X3 


ZR 

XOf GOTC 


Call AB6N0 


♦ IMTlALlZt PEGISTFRS FOR LOuPS 

% 


GCT»^ 


LPO^'iK 


lOQPJ 


SAl 

ARGLl 

A 

Sm2 

aRGL2 

0 

Sm 3 

ARG13 

C 

SOI 

-1 

-1 

5 BA 

PI 

-1 

SB7 

XA 

JJ 

SA4 

MAX 


s«? 

XA 

MAX 

5AA 

MIX 


$03 

XA 

MIX 

SAA 

IT 


$B*^ 

XA481 

iI-1 

SAA 

KK 


jXC 

XA 

KK 

SOA 

BA-Bl 

I»I41 

SX7 

Xl+BA 


SR6 

00 

J«0 

sxe 

BO 


SAA 

X7 

A 

5 AS 

X2^86 

P 

FXA 

XA^XS 

A*P 

Fxe 

X6 + XA 

♦ C 

5X7 

X7 + 82 

STEP A 

SB6 

B6-8i 

3 « J4>1 

LT 

P6f «7f LOQPJ 


NX6 

X6 


SA6 

X3 + BA 

-c 

LI 

PA,B5,100PIK 


SOA 

01 

I-O 

sxo 

XO + Bl 

K*K + 1 

wX2 

X2 + 83 

STEP B 

SX 3 

X3 + B2 

STEP C 

NZ 

XOf LOQPIK 


EG RETURN 
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lOFNT SUf-ULT 

* 

CLhMENT HAlftlX KULT C-(A*)*9 ♦ C LCWfcf; TPIAK'GILAP 

* 1 USQO kONS PACKWAPDS PROH Ii TQ 1 TO SAVE A B PEC-ISTFR 

* 



u:>E 

DATA 

HA X 

OSS 

IB 

KK 

BSS 

IB 

n 

«ss 

IB 

HI X 

Bw5 

IB 

I 

BSS 

18 

J 

BSS 

18 

BLCc 

8SS 

IB 


FNTRYP SOhULT 


♦ 


♦ SAVE 

EOPHAl 

PAPAHETPRS 


SA2 

Al + IB 


Sa3 

A2 + 1R 


8X6 

XI 


8X7 

X2 


SA6 

ALOC 


SA7 

BLOC 


PX6 

X3 


SA6 

CLHC 


CALL 

GET# (0#HAX#KK#iI)#ALOC 


CALL 

GET#(0#MX#1#J)#CLCC 

♦ 




♦ initialize REGIiTERi FDP LOOP 

♦ 



SAl 

PAX 

X1*H.AX 


SA3 

HIX 



SAA 

KK 



SA5 

II 



SB7 

XA 

B7*KK 


TX6 

X3AX5 



SBA 

X6 

BA-HIXAII 


SB3 

X3 

R3-HIX 


SRI 

-1 

81—1 


SA3 

BLOC 



1X3 

X3-X1 

X3-ADDK(B)-MAX 


SB2 

X3 

R2«X3 


SA2 

CLDC 



SAA 

ALTC 



SX5 

xs+ei 

(iI-1) 


1X2 


X2-A00P (C(n#l)> 


IXC 

X5AX1 



TXO 

XO^XA 

XO«ADDR(Aa#n)) 


SB5 

0 


lOOPTJ 

SA5 

X2 + B5 



Bxe 

X5 



see 

BO 



SR2 

X1^B2 

STEP B COLUMN 

ldopk 

SAA 

X0 + B6 

A(K#I) 


SA5 

B24B6 

8(K#J) 


PX5 

XAAX5 



FX6 

X5 + X6 

♦C 


NX6 

X6 



SB6 

B6-B1 

K-K41 


LT 

B6#B7» LCOPK 



SA6 

X2^B5 

•CCI* Jl 


SB5 

B5 + B3 

STEP C COLUMN 


LT 

85#8A#L00PIJ 



SBA 

BA-83 

REDUCE J LIMIT 


SB5 

80 



IX^^ 

xo-xi 

STEP A COLUMN BACKWARDS 


SB? 

X3 



SX2 

X2-1 

STEP C RCW BACKWARDS 


LT 

BOf BAfLOOPIJ 


• 

EO return 



END 
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lOENT 7HULT 


COMMENT C-C TAKEi ADVAKTAGF CF STRUCTURE OF A 


USE DATA 


max 

BSS 

IB 

IT 

BSS 

IR 

JJ 

ess 

le 

MIX 

ess 

IR 

J 

BSS 

18 

KK 

BSS 

IB 

SAVEC 

BSS 

IB 

ONE 

DATA 

1.0 


ENTRYP ZMULT 

♦ 

♦ SAVF FORMAL PARAMETERS 


SA? 

A14>1B 

SA3 

A2 + IB 

BX6 

XI 

RX7 

X2 

SA6 

ARGLl 

SA7 

AR6L2 

BX6 

X3 

SA6 

SAVEC 


CALI GFT>(0^MAX»II>JJ)»AR6U 
CALL GET#(0»MIX,JfKK),ARGL2 
* COMPARE VALUES OF J AND JJ 
SA3 J 
SAA JJ 
IXO XA-X3 
ZR XO^GOTC 
CALL ABEND 


♦ INITIALIZE REGISTERS FOP LOOPS 


0 




GOTO 

SBl 

-1 

81—1 


SA2 

MAX 



SB2 

X2 

B2«MAX 


SAA 

KK 



TXA 

X2«XA 



SB7 

XA 

B7-KKRMAX 


SA5 

jJ 



SX5 

X5 + B1 



SA3 

ARGL2 



1X7 

X3^X5 

X7-L0C(B(JJf in 


1X5 

X2^X5 



SB5 

X5 

B5-(JJ-1)*MAX 


SA5 

iT 



SA2 

X54B1 

A2-II-1 


SB6 

B5 + A2 

B6*(JJ-1)^MAX 4 


SA3 

SAVEC 



S83 

A2 



SX2 

X3 + B3 

X2«L0C(C(1I»1) ) 


PXO 

X2 

X0-X2 


SAB 

MI X 



SB3 

X3 

B3-MIX 


SAl 

ARGLl 



SAl 

XI 

A1*L0C(A(1»1)) 


SAA 

ONE 



NXl 

XA 


LCGPTJ 

SBA 

60 

INITIALIZE IK 


SAB 

A1 + R6 

AIJ 


ZR 

X3»SA0 



SA5 

X7 

INITIALIZE B (JK • JJ) 


TXA 

X3-X1 



NZ 

XA#L0uP2 


LOOPl 

SAA 

X2+BA 

C(1K) FRCM MEMORY 


FX6 

XAfX5 

C 4 B 


NX6 

X6 



SA6 

X2+BA 

C TC MEMORY 


SBA 

6A4B2 

1K»IK4MAX 


SA5 

A5 + B3 

JK-JK4MIX AND B(JK) FRCM MEMORY 


LT#I 

BA>B7, LOOPl 


SAO 

SX2 

X2481 

DECREhENl ROV OF C 


S86 

B640I 

DECREMENT ROW OF A 


GE 

B6>B5> LGOPIJ 



2.7 


♦ 


♦ 

LOOP2 


♦ 


S85 B5-82 DECREMENT COLUMN OF A 

SB6 A2^B5 

SX7 X7+B1 DECREMENT ROW OF B 

BX2 XO RESET ROW OF C 

GE B3,B0#IQ0R1J 

EQ RETURN 

SA^ X2+BA CIIK) FROM MEMORY 

FX6 X3*X5 A*B 

FX6 X64X4 +C 

NX6 X6 

SA6 X2+BA C TO MEMORY 

SBA B4+82 IK-IK+MAX 

SA5 A54B3 JK.JK+HIX AND B(JK) FROM MEMORY 

LT>8A#B7»L00P2 

EO SAO 


END 



IDENT MACROS 
STEXT 

♦ 

ENTRYP 

MACRO NAME 
ENTRY NAME 
USE 0 

TRACE. 

WFD A2/7L.NAME 
VFO IB/NAME 

SAVEAO 

BSS 1 

RETURN 

SAB SAVEAO 
SAD X5 

NAME 

EO ♦♦400000B 
SX6 AO 
SA6 SAVEAO 

ENTRYP. 

EOU NAME 

NCALLS. 

SET 0 
ENOM 

A 

CALL 

MACRO SU6>ARGL»ARGNAM 
LOCAL NONAM 

NAMED 

IFC NE»AARGNAMA^ 
SAl ARGNAM 

ARGLOEF 

IF -DEE^ARGNAM 
USE ARGLIST5 

ARGNAM 

BSS 0 

NAMED 

ELSE 

ARGLOEF 

IFC NE,^AR6L^* 
SAl NONAM 
USE ARGLISTS 

NONAM 

BSS 0 

NAMED 

ENDIF 
IRP ARGL 
VFO 60/ARGL 
IRP 

BSS7 1 
USF • 

ARGLOEF 

ENDIF 

NCALLo. 

SET NCALLS. +1 

♦ 

PJ -X.SUB 

“ 

VFO 12/NCALLS.#16/TRACF 
FNOM 

A 


END 


23 
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2 . 8 Small Computer Systems 

This section discusses the implementation of the MMLE3 program on small 
computer systems . The program is oriented toward large systems with such soft- 
ware features as NAMELIST input and CalComp plotter routines . Scratch disk space 
is assumed to be abundant . Finally , the program structure is optimized for seg- 
mented or overlaid loading , without which the core requirements would be quite large . 

The essential heart of the program , however , is fairly compact . By sacrificing 
some versatility and convenience features , the program can be reduced to the point 
that it will run on a minicomputer , Several suggestions for such reductions are 
made below . The specific reductions chosen will depend on individual requirements 
and capabilities . Several of the suggestions involve eliminating subroutines . This 
can be accomplished either by actually removing the subroutine and all calls to it, 
or by replacing the subroutine with a dummy consisting just of a return. 

Several of the suggestions will be recognized as replacing overlaid loading by 
a less automated equivalent— dividing the program into separate jobs run in the 
same core. 


2.8.1 Matrix Dimensions 

The most obvious saving is in reducing the maximum matrix dimensions as 
discussed in section 3.2. The dimensions in the published program are larger 
than needed for many applications , particularly applications appropriate to mini- 
computers . 


2.8.2 Predicted-Derivative Input 

The predicted-derivative input routine WTIN is the first candidate for reduction . 
In a nonoverlaid load , WTIN is extremely wasteful of core . The simplest solution 
is to eliminate WTIN from the MMLE3 program. WTIN can be run (if it is needed at 
all) as a completely separate program simply by replacing the subroutine card with 
a program card and replacing the return with a stop or by writing a driver that does 
nothing but call WTIN . The only communication between WTIN and the rest of the 
program is the file UWT created by WTIN . This file can be saved on disk or copied 
to cards . The UWT file can also be easily created without using WTIN for many 
applications. Reference 1, section 4.2.2, describes this file. 

An alternative to eliminating WTIN is to reduce the size of the 20 by 20 matrix 
in WTIN and/ or store it in common blocks used elsewhere in the program for other 
purposes. 

For larger reductions, the UWT file can be eliminated. Subroutines WTIN, 
WTDEF , and WTTRAN would then be unneeded , as well as the buffers and disk space 
for the file . 
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2.8.3 


2.8.3 Plotting 

Time history plots are by far the most useful "frill" of the MMLE3 program . In- 
deed , this frill is almost a necessity for analyzing real flight data . The plotting 
routines and associated system software use a large amount of core in a nonoverlaid 
environment . At the cost of some inconvenience , the time history plots can be cre- 
ated by a separate job . 

To remove the time history plotting from the MMLE3 program , delete the sub- 
routine THPLOT and the call to PLOT at card MMLE3 .71 . Subroutines AXES , LINES , 
PLTDAT, SCALE2, SYMBL4, TITPLT, and the CalComp routines are all called from 
THPLOT , and are thus also eliminated . 

Subroutine THPLOT expects the time history data on the scratch file UT2 and 
information describing the plot options in common blocks COM , HEADING , INTEGR , 
SIZE , and TOPLOT . The values in blocks FILES and MAXIMS must also be defined . 

If THPLOT is run as a separate program , it is probably more convenient to redefine 
the options and information in the common blocks than to communicate it from MMLE3. 
The time history data of file UT2 can be stored on disk or tape . Alternately , these 
data could be punched on cards (in which case file UT2 is not needed in MMLE3) . 

If disk or tape storage is unavailable and the number of cards required impractical , 
the data of file UT2 can be recreated by a driver routine for THPLOT by using the 
final estimated system matrices . This would require a moderate amount of coding 
most easily done by cannibalizing subroutine GIRL , throwing out parts used only 
for the gradients . 


2.8.4 State Noise Option 

If the state noise option is not needed , a significant amount of core can be saved . 
The state noise option cannot be run as a separate job like the wind-tunnel input or 
time history plotting , so removing this option constitutes a sacrifice of capability , 
not just of convenience . The state noise option is removed by eliminating subroutines 
FADJ , FLIMIT , KALMAN , and GRADE . The routines ADDPAR , GETPAR , GRADP , 
LYAPCB, MOVE, MULTT, RICATC , SSIMEQ, and TRANSP are called only from the 
four primary state noise routines (directly or indirectly) and can thus also be 
eliminated. Common block GRDCOM and cards 71 to 82 should be removed from 
subroutine GRAD . The routines EIGENG , BALANC , HQR , HQR2 , ELMHES , ELTRAN , 
and BALBAK can be eliminated if the state noise option is removed and the computation 
of the determinant at cards RESIDE . 56 to 62 is eliminated . Several matrices used 
only with state noise (F , FV , APRF , P , and KGAIN) could be eliminated along with 
short sections of code in the rest of the program , but such savings are not large . 


2.8.5 Miscellaneous Files 

Several files used by the program can be eliminated to lower disk utilization . 
Elimination of files UWT and UT2 is discussed in sections 2.8.2 and 2.8.3, respectively . 
File UTHOUT is an easily sacrificed frill written by routine THOUT . The punch file 
UPUNCH can also be dispensed with, in which case subroutine OUTPUN would be 
eliminated . 
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2 . 8.6 


The file UDATA , read by subroutine READTH , is not used if the time history 
is read from cards . Alternately , the input time history file can be prepared before- 
hand , with times selected out and data corrections made , and attached directly as 
file UTl. In this case most of subroutine THDATA can be eliminated; only the 
averaging and the call to user routine AVERAG need remain . If absolutely no disk 
or tape storage is available , the time histories can be stored in core , but this uses 
a large amount of core even if one is very selective about the number of signals used . 
If there is a capability for having only single file on disk or tape , that file should 
be UTl. 


2.8.6 HQR and HQR2 

Subroutine HQR is the same as HQR2 without the eigenvector computation . 
Therefore, HQR can be eliminated by adding an argument to HQR2 to control skip- 
ping the eigenvector computations . The change is easy and saves a worthwhile 
amount of core if segmentation is not used. The two routines are separated in MMLE3 
because it was desirable to use the EISPACK routines as published in reference 6 , 
and because there is no penalty for the code duplication if segmented loading is used . 

If the state noise option is removed , the routines HQR and HQR2 will normally 
also be removed and this question will be moot . 


2.8.7 Minimum Program 

If MMLE3 is to be run on a small minicomputer , it may be necessary to make 
further reductions than those described above . This section describes the reduction 
of MMLE3 to about the minimum usable size. 

The actual derivative estimation is performed by subroutine NEWTON and the 
subroutines called by NEWTON . The rest of the program can therefore be dispensed 
with , except for the minimum needed to define quantities used by NEWTON . In 
particular, subroutines WTIN, TITLES, EDIT, MATSET, MTLOAD, THDATA, MATDEF, 
COMP AT , ALLOW , SUMOUT , THPLOT , and subroutines called only through these 
routines can be eliminated. Section 4 describes the general relationships of these 
routines . Figure 1 can also be helpful in checking which routines and common blocks 
can be eliminated in this process . Some type of minimal input routine will be needed 
to define the common block parameters used by NEWTON . In addition to the basic 
matrices and options , the data in common block DETERM must be defined . The 
unmodified MMLE3 program reads the data in a more convenient format , and then 
subroutine ALLOW creates the lists in block DETERM based on the input . This 
makes the program easier to use , but uses significantly more memory than reading 
the data directly in the format required by NEWTON . 

A minimal size MMLE3 would also incorporate the previously discussed modi- 
fication to eliminate the state noise option to reduce the size of the part of the 
program under subroutine NEWTON . A few more words could also be saved by 
eliminating the separate storage of the M suffix matrices in common block MATRAT . 
The M suffix matrices could be equivalenced to the matrices in block DIMMAT 
instead of stored separately . 
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3.0 


3 . 0 MATRIX STORAGE 


This section discusses the matrix storage conventions used by the MMLE3 
program . Methods of changing maximum matrix dimensions are also shown , 


3 . 1 Conventions 


The MMLE3 program must be able to work with matrices of varying size without 
recompilation . The FORTRAN language is not well equipped to handle variable size 
matrices , A set of matrix storage conventions has been adopted to partially alleviate 
this problem , Figure 2 illustrates the conventions discussed in this section . 


Number of physical columns 




Number of 
physical rows 


Matrix locations used 
(logical size) 




Unused 





Flag 

Number of 
rows 

Number of 
columns 

Name 

Unused 



Figure 2. Matrix storage conventions. 


27 



3.1 


The numbers in the FORTRAN dimension statement are referred to in this report 
as the physical dimensions of a matrix. These values control the amount of space 
allocated for the matrix by the FORTRAN compiler , and thus the largest usable di- 
mensions of the matrix . Changes of the physical dimensions of a matrix require 
recompilation and are discussed in section 3.2. 

Although the physical dimensions of a matrix are fixed , the program can store 
a matrix smaller than the physical dimensions as a partition of the physical matrix, 
ignoring the remaining locations . The size of this partition in use , referred to as 
the logical size of the matrix , may vary during program execution . No special 
conventions have been adopted for vectors (unless they are stored as a matrix with 
one logical column or row) . The following conventions are used in MMLE3 to keep 
track of the physical and logical dimensions of matrices . 

The last physical row of each matrix is reserved for information about the 
matrix; therefore , the physical number of rows of a matrix must always be at least 
one more than the logical number of rows . The physical number of columns of 
every matrix should be at least 4 , because 4 locations in the last physical row 
are used. 

The first location in the last physical row is used as a flag to find the row . 

The value used as a flag is defined in subroutine GETSET . Currently this value 
is the real number equivalent to the characters "TEST" left -justified and blank-filled. 

92 

(On a CDC 60-bit computer, this value is approximately 3.149 X 10 .) The flag 

value should be one that will not occur in any of the matrices stored. Note that the 
chances of any randomly picked value occurring in any given location are about 1 
9 18 

in 4 X 10 on a 32-bit computer and 1 in 10 on a 60-bit computer . The odds in 
practice are even more unlikely for very large values such as the one used here , 
because the units used for most physical problems tend to keep values within a 
couple of orders of magnitude of 1 . On a CDC computer , the value "negative indef- 
inite" is a particularly good choice , because it csin never be the result of an 
arithmetic operation . (Subroutines GET and GETP would have to be slightly modified 
if "negative indefinite" were used.) 

The second and third locations in the last physical row are the number of logical 
rows and columns , respectively . These values are stored as real numbers . The 
fourth location in the last physical row contains the name of the matrix in Hollerith 
format . This name is used when the matrix is printed or punched . The name DONT 
is reserved for a special convention . The matrix print or punch routines (SPIT 
and PLOP) will not output matrices with a name of DONT . This convention simplifies 
control of which matrices will be output for a given case . 

The functions of all the matrix manipulation routines are described in appendix A , 
but a few of the most important will be briefly described here . GETSET is an initial- 
ization routine that must be called before any of the other matrix routines . Routines 
SET , SETl, and SET2 set elements in the last physical row of a matrix. GET and 
GETP retrieve information about the physical and logical size of a matrix. UNSET de- 
letes the flag from the first element of the last physical column of a matrix (this allows 
the same space to be reused as if it were a matrix with different physical dimensions) . 
ABEND is an error routine; it intentionally causes an end-of-file error on the card 
reader file in order to get an error traceback . 
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3.2 


3.2 Changing Maximum Dimensions 

The matrix storage conventions allow the logical dimensions of a matrix to be 
changed within limits without recompilation . However , if a system to be analyzed 
exceeds the physical dimensions of a matrix, the program must be recompiled with 
larger dimensions. Conversely, if most of the analysis at an installation is done on 
very small systems , it may be prudent to recompile the program with smaller dimen- 
sions to avoid wasting core . 

Changing the physical dimensions in a program as large as MMLE3 can be a 
major task. The coding of the MMLE3 program was specifically designed to simplify 
this task. Every card that must be changed in order to change the physical dimen- 
sions is in a common deck (see sec. 1) . This convention minimizes the number of 
places to check for changes . 

In order to simplify the task further , program COMPUN was written to punch 
out a complete set of the common decks with altered dimensions . The only input 
to COMPUN is a NAMELIST defining the desired matrix sizes. A listing of this 
program is shown in appendix C . The punched output from COMPUN can be used in 
the COMSUB program (appendix B) to create the actual FORTRAN code. Alternately, 
CDC UPDATE (ref. 2) can be used in place of COMSUB . Users of UPDATE will find 
it convenient to change COMPUN so that it punches out a correction set . This is 
done by changing each *COMDECK card in the data for COMPUN to an appropriate 
delete card; an ident card will also be needed at the beginning, and resequence 
cards may be desirable at the end (ref. 2) . 


4.0 PROGRAM STRUCTURE 


This section gives a brief overview of the structure of the MMLE 3 program . 

Because of the large number of subroutines, it is difficult to understand the oper- 
ation of MMLE 3 by studying the individual subroutines. A guide to how the sub- 
routines fit together as a coherent whole is needed. With this guide, the user should 
be able to determine which subroutines are involved in any particular task; appendix A 
and supplement 1 can then be consulted to find details of how the task is accomplished . 

Figure 3 constitutes a structural diagram of the major routines of MMLE 3. 

Each box represents a task performed by one subroutine or a closely related group 
of subroutines . The primary subroutine name is underlined at the top of the box , 
and a brief description of the task follows . Important secondary subroutines called 
in the performance of each task are listed at the bottom of the box. Minor subroutines 
such as those for matrix manipulation are not included in the figure . The matrix 
subroutines are called throughout the program . 
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APPENDIX A 

DESCRIPTIONS OF SUBROUTINES AND COMMON DECKS 


This appendix contains descriptions of the function and operation of each of the 
subroutines of the MMLE3 program. The variables in each of the common decks are 
also described . The user routines described are the standard aircraft routines 
(see ref. 1 , sec, 4) . 

Supplement 1 contains microfiche listings of the subroutines and common decks; 
it should be consulted in conjunction with the descriptions in this appendix . The 
common decks are listed first. Next, all of the *CALL cards are listed separately. 
The separation of these cards from the FORTRAN listings is inconvenient , as cross- 
checking between the list of the *CALL cards and the FORTRAN listing is sometimes 
necessary , but a more convenient format of the listings is not easily obtained . The 
following pages contain lists of the UPDATE correction idents and deck names , com- 
pleting the listings produced by UPDATE . The largest portion of supplement 1 
consists of the FORTRAN listings and reference maps . The EISPACK routines are 
listed separately at the end of supplement 1 , preceded by an UPDATE deck list for 
these routines . 


A . 1 Common Decks 

The common decks will be described in the order that they are shown on the 
listing of supplement 1 . The common decks of the basic program are described 
first, followed by those of the standard aircraft routines. Following each description 
is a listing of the common deck . 


A , 1 , 1 Basic Program 

The first common deck of the basic program is common deck HISTORY; the 
remaining common decks are in alphabetical order . 

A . 1 . 1 . 1 HISTORY . - Common deck HISTORY contains only comment cards , This 
common deck is intended to give a history of program modifications . For each modi- 
fication to the program , one card briefly identifying the modification should be added 
to this common deck. Common deck HISTORY is called only from the main program. 


♦C^MDECK HISTORY 
C MODlFICAtlUH HTSTDPY: 

* END OF MODTFICATnNS. 


HISTORY 1 
HISTORY 2 
HISTORY 3 
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A. 1.1. 2 


A . 1 . 1 .2 AMCOM . - Common deck AMCOM contains the vector BCONST . This 
vector is defined by subroutine CALLAM; it contains the dimensionalization ratios 
from the M suffix matrices . For dependent variables in hard constraints , the corre- 
sponding element of BCONST is the constraint ratio times the dimensionalization 
ratio from the M suffix matrix . 


*C^MnECK ihCOf 

CJM1C1M /AMCOH/ I^CINSK 50) 


AHCOX 1 

AftCOH Z 


A. 1.1. 3 AVGCOM. - Common deck AVGCOM contains the average values of the 
measured states (ZAVG) , controls (UAVG) , and extra signals (EXAVG) . The 
averages are computed for the entire , dimensioned lengths of these vectors , regard- 
less of the vector size used in the system equations . The averages are taken over 
all of the maneuvers if multiple maneuvers are used . Standard deviations , maxima , 
and minima are also stored in AVGCOM . The relative position of the variables in this 
common block should not be disturbed, since subroutines THDATA and DlMl depend on 
this order, treating ZAVG, UAVG, and EXAVG as a single, concatenated vector. 


‘■“□IDECK AVGCOM AVGCOM 1 

COMMON /AVGCOM/ /AVGt 0B),UAVG( 0A),EXAVG( 20)»?SI6( 08), AVGCOM ? 

USIG( 0A),EVjI6( 20),ZHiMM( 08),UM1NM( 0«),EVMINM( 20), AVGCOM 3 

7MAXM( O0),'JMA*M( uA),EXMAVM( 20) AVGCOM A 


A . 1 . 1 .4 BIASES . - Common deck BIASES contains the bias vector UOFF and YOFF 
defined by subroutine INIT . The role of these biases is discussed in reference 1 , 
section 3.1. 


• •“OMDECK BIASES 

COMMON /BIASES/ UOFFI ,;A),Y0FF( 08) 


BIASES 1 

BIASES 2 


A. 1.1. 5 BILIN . - Common deck BILIN contains the measured time history data 
for one time point. The observations (Z), controls (U) , and extra signals (EXTRA) 
are obtained from the time history data file. The bias vector (ONES) is computed 
in subroutine GIRL (lines 61, 62, and 67) as described in reference 1, section 3.1. 
The logical variable TIMVAR is described in reference 1, section 3.3.8(20) . Average 
values from common block AVGCOM are placed in Z , U , and EXTRA by subroutine 
DIMl when USEAVG is TRUE . Subroutine GIRL controls USEAVG. These average 
values are used in the computation of the Kalman gain matrix and its gradients . The 
relative position of Z , U , and EXTRA in the common block should not be disturbed , 
or the average values will not be placed correctly . 


♦COMDECK BILIN BILIN 1 

COMMON /BILIN/ USEA V6,TIMVAR, Z ( 08), U( OA), EXTRA! 20), ONES! CA) BILIN 2 

logical USEAVG, TIMVAR BILIN 3 
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A. 1.1. 6 


A . 1 . 1 . 6 COM . - Common deck COM contains information about the time points 
used for a case . NCASE is the number of maneuvers being analyzed (ref. 1, 
sec. 3.3.8(1)) . NPTS contains the number of time points in each maneuver, and 
NPTT is the total number of time points for all of the maneuvers . ITMSTS contains 
the total time in milliseconds of the first time point in each maneuver , 


‘CIMOECK COH 

cnimoh /COM/ NCASE»NMTT,NPTS(15»,ITMSTS(15) 


CON 1 

COM ? 


A. 1.1.7 CRMAT . - Common deck CRMAT contains the Cramer-Rao bounds for the 
nondimensional matrices . 


♦COMPECK C»NAT CRMAT 1 

COMMON /CRMAT/ AC ( 08, 07),aC( Ctl, 0A),SC( 08, 0A»,RC( 08, 07), CRMAT 2 

CCl 07, 07),0C< 09, 0A),HC( c9, 0A),EC< 09, 07),EC( 08, 07) CRMAT 3 


A . 1 . 1 . 8 DETERM . - Common deck DETERM contains information describing the 
unknowns to be determined. The first half of subroutine ALLOW defines the vari- 
ables in this common deck . NVAR is the total number of unknown locations . It 
equals the number of independent unknowns plus the number of constraints; thus , 
some locations may be counted more than once if they are dependent variables in 
more than one constraint . The remaining variables in the common block are vectors 
of length NVAR. IMAT contains the matrix number. The matrices (ref. 1, 
sec. 3.3.11(1)) corresponding to the numbers in IMAT are 1 = AN, 2 = BN, 3 = SN, 

4 = RN , 5 = CN , 6 = DN , 7 = HN , 8 = EN , 9 = FN , and 10 = initial condition . IROW and 
ICOL contain the row and column numbers , respectively (for unknown initial condi- 
tions , the column number is ignored) . ILOC gives the location of the associated 
gradient variable, and ACONST is the constraint ratio. Each independent variable 
will be associated with a unique gradient variable and the constraint ratio will be 
defined as 1 . Dependent variables in hard constraints will be associated with the 
same gradient variable as the corresponding independent variable; the appropriate 
constraint ratio will be used . The information in common block DETERM is printed 
in the output under the heading "LOCATION INDICES ." 


♦CTMPFCK OFTERM PETERM 1 

CIMMQN /DETERH/ NVAR,IMaT( ! 0 ),IR 0 R( 50),ICCL( 50),ILaC( 50), OETERI* 2 

AC0NST( 50) DETERM 3 


A . 1 . 1 . 9 DIMMAT . - Common deck DIMMAT contains the dimensional system 
matrices . They are defined in subroutine DIMl . 


♦'•PMDECK DIMMAT DIMMAT 1 

COMMON /DIMMAT/ A ( 08, 07), B ( 06, 0A),S ( 08, 0A),R ( 08, 07), DIMMAT 2 

C « 09, 07),P t 09, 09), H { 09, 09), E ( 09, 07) DIMMAT 3 
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A. 1.1. 10 


A . 1 . 1 . 10 DUMCOM . - Common deck DUMCOM contains three matrices used for 
scratch storage in several subroutines . 


*:nnnECK DU'^CCM 

CLillUN /OU.ICJM/ 0UM( Cbi C7lf0U»'2( 08» 07),DUH3( 08, C7» 


OUNCOM 1 

OUMCOT* 2 


A. 1.1.11 DUMVEC . - Common deck DUMVEC contains the vectors DUMX and 
DUM2 used for scratch storage. 


♦'•OMLECK DUMVEC Ol'MVEC 1 

CJMMuN /DUMVEC/ DUMX( 07),DUM2X( 07) POMVEC 2 

A .1.1. 12 ECOM . - Common deck ECOM contains matrices of the state equation 
multiplied through by i? ^ . The specific variables are RI = R ^ , RIA = R , 

RIB = R~^B , and RIS = R ^S. These matrices are defined by subroutine DIM2. 


♦CTMDcCK FCDM 

COMHON /ECOM/ RIA( 08, 07),RI6( OP, 04).RIS( 08, 04).»T( 08, C7) 


ECOM 1 

ECOM 2 


A. 1.1. 13 ERLIST . - Common deck ERLIST contains information on the convergence 
of the cost functional . NITER is the current interation number , and ERRVEC is a 
vector containing all of the previous values of the cost functional. The logical vari- 
able BLOWUP is set if the error becomes unreasonably large , indicating probable 
divergence . This causes iteration to stop prematurely . The dimension of ERRVEC 
can be changed without affecting any other program dimensions; this dimension 
limits the allowable number of iterations. 


♦'tlriDFoK PRLIST E«LIST 1 

ZmnDH /ERLIiT/ BLOWUP, MITER# ERRVBC(50) ERLIST ? 

LJGICtL 9LQWU0 ERLIST ^ 


A . 1 . 1. 14 FCOM . - Common deck ECOM contains the state noise power spectral 
density matrix , F . 


♦DTMOECK FCOM 

CnMMDN /FCO"/ F( D8, 37) 


FCOM 1 

FCOM 2 


A. 1 . 1.15 FILES ■ - Common deck FILES contains the variables used for I/O unit 
numbers. These variables are discussed in reference 1, section 3.2. Their values 
are defined by subroutine VARDEF . If these values are changed , the program card 
(cards MMLE3.2 to 6) must also be changed. 


♦OOMDECK FILES FILES X 

COMMLn /FILES/ UC 480 , UPUMCH ,U PR I NT.UDAU, UTX,UT2,U1MQUT,UWT, UPL OT FILES 2 

INTEGER UCARD,LIPUNCN,UPKlNT,unATA,UTX,UT2,UTH0UT,UWT,UPL0T FILES 3 
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A. 1.1. 16 


A. 1.1. 16 GICOM . - Common deck GICOM contains the GGI matrix and related 
information . RSQ and FRSQ are the sample covariances of the raw and filtered 
residuals (fit errors) , respectively . WRSQ and WFRSQ are the diagonal weighted 
fit errors. FREQCR, ITG, RLXG, and DIAGG are the input variables (ref. 1, 
sec. 3.3.8(23) to (25)) that control the residual filter and the G determination. If 
the residuals are not filtered (FREQCR equals 0) , FRSQ will be the same as RSQ . 
ERRFLT is the total filtered error sum , and SGNLS is the number of weighted sig- 
nals. FCl and FC2 are the constants computed to implement the residual filter. 


♦ 'TMOECK GICOM GICCM 1 

C1MM0N /GICOM/ ITG. DIAGG. FREOCR»RLXG. FCl. FC2»ERRFLT«SGHI.S. GICOM 2 

G6I( 09. 08).RS0( J9. 03).FRS0< 09. 08).RR50< 08).WFRS0( OS) GICPM 3 

LOGICAL OIAGG GICOM A 


A. 1.1. 17 GRADS . - Common deck GRADS contains the gradients of the states 
and observations . GRADX is the gradient of the state and GRADY is the gradient of 
the observation. The residual is stored as an augmented column of GRADY. GRADl 
is a scratch matrix used in the computation of both gradients . The matrices are 
also used for scratch storage in subroutine FLIMIT . 


‘COMPECf GRADS 

CTMMON /GRaOS/ GRA0X( 0f>. 35).GRADY( 09. 35).GRAD1( 09. 35) 


GRADS 1 

GRADS ? 


A . 1 . 1 . 18 GRAD$ . - Common deck GRAD$ contains vectors needed only in 
subroutine GRAD . XT 12 is the average of the state at the beginning and end of a 
sample interval. XDT2 is the derivative of the state at the end of the interval; 
XDT12 is the average of the state derivatives at the beginning and end of the sample 
interval. The XDT's are computed ignoring state noise. 


♦'•OdCECK GRA0$ 

OIMcnsipn XDT2( 07).XOT12( 07).XU2( 07) 


GRAOS 1 

GRAOS 2 


A . 1 . 1 . 19 GRDCOM . - Common deck GRDCOM contains the triply dimensioned 
array DK . The third index of the array corresponds to the list of unknowns that 
affect the K matrix . For each of these unknowns , DK ( • , • , i) contains the gradient 
of K with respect to that unknown . DK is defined in subroutine GRADK . 


♦CQ/IOECK GRDCOM 

COMMON /GRDCOM/ OKI OP. OP. 15) 


GRDCOM 1 

GRDCOM 2 


A. 1.1.20 GRSIZE . - Common deck GRSIZE contains information about the size 
of the gradient vectors. JKMMl is the number of independent unknowns. NK is the 
number of independent unknowns affecting the K matrix . 


‘COH05CK GRiI7E 

CjMMLN /GRSI2F/ JKhMl.NK 


GRSIZE 1 

GRSIZE 2 
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A. 1.1. 21 


A. 1.1.21 HEADNG . - Common deck HEADNG contains labels and titles. SIGLAB , 
XLAB , CONLAB , and EXLAB contain the labels for the observations, states, con- 
trols , and extra signals , respectively . Two words are allowed for each label . 

TITLE contains the title card for the case . ABATE and ATIME contain the date and 
time if available to the program . The relative position of SIGLAB , XLAB , CONLAB , 
and EXLAB should not be changed , as subroutine THPLOT depends on this relation- 
ship (card THPLOT . 157) . 


‘“JMljcCK HEADNG HEADNG \ 

GQMMDN /HEADNG/ TTTL E ( ?0 ), AOAT F , ATIME, HEADNG 2 

SiGLA6(E, 0e),XLfc8(2, 07 ) , CnNLA8( 2, 0A),EXLAR(2. 20) HEADNG 3 


A . 1 . 1 .22 ICOND . - Common deck ICOND contains information about the initial 
conditions. USERIC and VARIC are input variables discussed in reference 1, 
section 3.3.8(26) and (27) . VARICS is the Boolean sum of the elements of VARIC . 
DXIC is the initial condition increment estimated when elements of VARIC are TRUE . 


‘'•TMDFCK ItOND 

COMMON /ICDND/ UCEMC , VaR iCS, VARIC ( 07 ) , DX IC ( 07) 
LUGUAL U3ERIC,VARICi, VARIC 


IC3ND 1 
ICON!) 2 
ICOND 3 


A . 1 . 1 .23 INMAT . - Common deck INMAT contains information from a matrix 
header card during matrix input . ALAB is the matrix name , and IM is the corre- 
sponding matrix number defined by function MATNO . II and JJ are the numbers of 
rows and columns , respectively . 


•comdeck innat 

COMMON /INMAT/ ALAB,n,JJ,TM 


INMAT 1 

INMAT 2 


A . 1 . 1 .24 INOPT . - Common deck INOPT contains the logical variables CARD and 
TAPE, described in reference 1, section 3. 3. 8 (2). 


♦ C'lMDECK INQPT INOPT 1 

COMMON /INJPT/ CARD, TAPE INOPT 2 

LOGICAL CARO, TAPE INOPT 3 


A. 1 . 1 .25 INQRD . - Common deck INORD contains information about the order of 
signals on the time history data file . All of the variables are described in reference 1 , 
section 3.3.8(G) and (7). The relative position of ZCHAN , UCHAN, and EXCHAN 
should not be changed, since subroutine THDATA depends on this relationship. 


♦COMUECM INORD INOP.D 1 

COMMON /INORD/ NREC,ZCHAN( 0P),UCHaN( 0A),EXCHAN( 20) INORD 2 

Integer zchan,ughan,pxchan inord 3 
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A. 1.1. 26 


A . 1 . 1 . 26 INTEGR . - Common deck INTEGR contains data required for the integra- 
tion routine EAT . DT is the sample interval of the data (after any thinning) . NEAT 
is an input variable described in reference 1, section 3.3.8(15) . 


‘'"TIDECK INTEGR INTFGP 1 

raMNON /INTEGR/ OT»NEAT INTEGR 2 


A. 1.1, 27 KCOM . - Common deck KCOM contains the Kalman gain matrix , KGAIN , 
and the Riccati covariance matrix , P . 


♦ '“TNOECK KCON KCOR 1 

CCJM'iCN /KI.JN/ Rt 38» 07).KG*1M OE. 08) KCON 2 

REAL KGAiN KCON 3 


A. 1.1, 28 MAPCOM . - Common deck MAPCOM contains internal location maps . 
These maps are created by the second half of subroutine ALLOW and are used only 
in the state noise algorithm. Each vector in MAPCOM maps from a source list to a 
destination list. Each map vector in MAPCOM is the same length as the source list; 
each element in the map corresponds to an element in the source list . The value of 
each map element indicates the position in the destination list with which the cor- 
responding source element is associated . A value of 0 for any map element indicates 
that the corresponding source element is not associated with any element of the 
destination list . 

The variable names of the maps are all five characters , the first three of which 
are "MAP . " The fourth and fifth characters indicate the source and destination 
lists , respectively . The letters used for the fourth and fifth characters are U , G , 
and K . U represents the complete list of unknowns , including independent unknowns 
plus constraints; the length of this list is NVAR (common block DETERM) . G 
represents the complete list of gradient elements; the length of this list is JKMMl 
(common block GRSIZE) . Each independent unknown will correspond to one gra- 
dient element. K represents the list of gradient elements which affect the K matrix; 
the length of this list is NK (common block GRSIZE) . 


♦CTiOECK NaRClN 

COMiOS /NARCnN/ NAPUK( 50)»MARK6( 15) 


NARCO" 1 

)(APC0N 2 


A. 1.1. 29 MATIN . - Common deck MATIN contains matrices that are only used in 
the input section of the program. The information from these matrices is put into 
other forms for use later in the program . The matrices in this common block are 
the V suffix, APR prefix, and hard constraint (HARD) matrices. All of these matrices 
are described in reference 1, section 3.3.11(3) , (4) , and (6) , respectively. 


♦ '•nilLECK NATIN NATIN 

COHNON /NATIN/ AV( OP, 07),BV( 08, 0A),SV( 08, 0A),RV( 08, 07), MATIN 

CV( 09, 07),0V( 09, uA),HV< g9, 0A),EV( 09, 07),FV( 08, C7), NATIN 

APRA( 0 », 07),APRB< 08, 0A),APR5< 08, 0A),APRP( 08, 07), NATIN 

APRC< 09, 07),APRn( 09, 0A),APRH« 09, 0A),APRE( 09, 07), NATIN 

APRF( OP, 07),HAR0< 36,7) NATIN 


1 

2 

3 

4 

5 

6 
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A. 1.1. 30 


A. 1.1 .30 MATRAT. - Common deck MATRAT contains the M suffix matrices . 


♦ nMDPCK MATRAT 1 

CJMHON /KATRAT/ AM( Ot. 07).RM( 08. 0A».SM( 08. 0A).RM( 08. 071. MATRAT 7 

CM( O'?. 07),DM( O'?. C<t).HM( O'?. OA).FM( 09, 07) MATRAT 3 


A. 1 .1.31 MATLAB . - Common deck MATLAB contains the list of matrix labels 
and read flags . NMATS is the length of these lists (currently 31) . LAB is the list 
of matrix labels . INFLAG is the list of read flags . Each element of INFLAG is 1 if 
the corresponding matrix has been read from cards; the element is 0 if the matrix 
has not been read . INFLAG is initialized to 0 , and NMATS and LAB are defined in 
subroutine MATSET . INFLAG is then altered by cards MTLOAD . 25 and CONIN . 28 . 


♦riiDECK MATLAB f'ATtAR 1 

CJMMQM /M4UAP/ NMAT^U40(31)» XNElAGOI) P*ATLAB ? 

real IaB hATLAP 3 


A , 1 . 1 . 32 MATRIX . - Common deck MATRIX contains the N suffix matrices , 
described in reference 1, section 3.3.11(1). 


♦''11JECK HATPiV matrix 1 

CQMflON /JIATkIX/ Aw( 38, 07),PN( u6, oA),SN( 03, 0A)>RN( 08, 07), f^ATRIX 2 

CN( )9, 07),0K( 09, CA),HN( 09, j^),EN{ 09, 07) MATRIX 3 


A . 1 . 1 .33 MAXCON . - Common deck MAXCON contains the maximum dimensions 
MAXHRD and MAXSFT for the constraint matrices, HARD and SOFT (ref. 1, 
sec. 3.3.11(6) and (7) , respectively) . The variables MAXHRD and MAXSFT are 
defined by common deck VARDEF$. Maximum matrix dimensions are discussed in 
section 3 and in reference 1 , section 2 . 


♦ tiAXCON 

C /MAXCON/ MAXHR0,MAXiFT 


MAXCON 1 

MAXCON 2 


A. 1.1. 34 MAXIM . - Common deck MAXIM contains some maximum matrix dimen- 
sions. NI , MAXTV, and MAXKV are defined by common deck VARDEF$. MAXXl and 
MAXZl are defined as MAXX + 1 and MAXZ + 1 in subroutine VARDEF. Maximum 
matrix dimensions are discussed in section 3 and in reference 1 , section 2 . 


♦CI'^ineCK MAXIM 

COiMON /MaXIM/ HAXX1,«AXZ1,NI,»<AXTV, MAXKV 


MAXIM 1 

MAXIM 2 


A. 1 .1.35 MAXIMS . - Common deck MAXIMS contains some maximum matrix 
dimensions . MAXX , MAXZ , MAXU , MAXB , and LEX are defined by common deck 
VARDEF$. LORD is defined as MAXZ + MAXU + LEX in subroutine VARDEF. Maximum 
matrix dimensions are discussed in section 3 and in reference 1 , section 2 . 


*''l' 10 ECK MAXIMS 

CIMM^N /maxims/ MAXX,hAXZ,MAXl,«AXP, LEX , lord 
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MAXIMS 

MAXIMS 


1 

2 



A. 1.1. 36 


A . 1 . 1 . 36 MODCOM . - Common deck MODCOM contains the logical variable UMOD . 
UMOD is TRUE if user routines are used; otherwise , it is FALSE . UMOD is defined 
by the main program; the input to control this is described in reference 1, sec- 
tion 3.3.2. 


♦■1.10ECK 100C0M KCDCPM 1 

CIMHGN /MODCOM/ UMOn MpocDM ? 

L0GiC4l UM1D MODCOM 3 

A >1.1. 37 OBSRV . - Common deck OBSRV contains matrices used for the computed 
observations and their gradients: ERIAC = ER~^A + C, ERIBD = ER~^B + D, 

ERISH = ER + H , ERI = ER ^ . These matrices are defined by subroutine DIM2 . 


♦ OBSRV 1 

/O&SRV/ t*?UC( 07)fERlR0( 09» 0^)»€RISH( 09, 0^), OBSRV ? 

09, 07) OBSRV 3 


A. 1.1. 38 OUTOPT . - Common deck OUTOPT contains variables controlling out- 
put options. All of these variables are described in reference 1, section 3.3.8. 


♦CTMdeCK uUTuPT 


/JUI JMT/ HKUtTX, ^RlMTYf ^RIhT3,PLaTEH,PUNCH#TPSTtPLTMAX# 
ERRTH 

logical PRINTX, PPTnT Y, print 0, PICT EM, PUNCM, TEST,EPRTH 


OUTOPT 1 
OUTCPT ? 
OUTOPT 3 
CUTCPT A 


A. 1.1.39 PBCOM . - Common deck PBCOM contains PB, the vector of changes in 
the coefficient estimates . The call to MVMULT at card NEWTON . 66 defines PB . 

It may be modified by subroutine FLIMIT . 


♦'TMPECK iJBCnf 

COMMON /BBCOM/ pp( 35) 


PBCOM 1 

PBCOM ? 


A . 1 . 1 . 40 PHICOM . - Common deck PHICOM contains the transition matrix and 
several products involving its integral . The transition matrix , PHI , is the expo- 
nential of R ^AAt. Call the integral of the transition matrix, xj/ . Then, 

PSIB = i|/R , PSIS = v|/R , and PSl = x(rR ^ . All of these matrices are defined by 
subroutine REAT . 


♦ '^QMuECK PHIcr^M 

common /phICCM/ pHK 08, C7),PSl( 06, 07),PSIB( 08, yA), 
PSI5( 08, 04) 


PHICOM I 
PHICUh 2 
PHICOM 3 
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A. 1.1.41 


A. 1.1. 41 RECRD. - Common deck RECRD contains one record from the time 
history data file . EOFTH can be set to TRUE by user routine READTH to indicate 
an end-of-file . T is the time in integer hours , minutes , seconds , and milliseconds . 
RECORD is the data for that time. The dimension of RECORD can be changed without 
affecting any other program dimensions . 


RECfcO 

CJflMGN /RECRO/ T ( ^ ) * RET Gfi C ( 100) 

LOGICAL PUPTH 
inEGtO T 


RECRO 1 

RECRO 2 

RECRO 3 

RECRO 4 


A. 1.1. 42 RICCOM . - Common deck RICCOM contains matrices used to compute the 
Kalman gain matrix and its derivatives. RIF = R , RIFRIF = R V(R V)*, 

CTG = C*(GG*)"^, and RIAP = R~^AP . DUMXZ and DUMZX are scratch storage. 


♦ flMOECK RICCOM RlCCn>> 1 

COMMON /RICCOM/ 0UM*7( 0R> v)<3, IT), RICC’M 2 

RiAf>( 08. 07),CTG( C6. ORI.RIFI 07I»RIFR1F( OF. 07) RICCOM 3 


A. 1.1. 43 SIZE. - Common deck SIZE contains the system vector sizes. MX is the 
length of the state vector , MZ the observation vector, MU the control vector, and 
MB the bias vector . These lengths are defined by subroutine COMP AT . They are 
discussed in reference 1, section 3.3.8(11) to (14). 


♦ CTMOtCK SIZE 

COMMON /Size/ MX.MZ.MU.M8 


SIZE 1 

SIZE 2 


A. 1.1. 44 SOFCOM. - Common deck SOFCOM contains the matrix of soft constraints, 
SOFT. 


♦01M0CCK «0FCGM 

COMMON /S1 = C0M/ iOFT( 11,7) 


SOFCCM 1 

SOFCOM 2 


A . 1 . 1 . 45 SUMCOM . - Common deck SUMCOM contains the SUM matrix , the second 
gradient of the cost functional . Only the lower triangular and diagonal parts of this 
symmetric matrix are stored . The first gradient is augmented as a last row or 
column. JKM is the logical dimension of the SUM matrix, i.e. , the length of the 
gradient vector JKMMl (common block GRSIZE) plus 1. In subroutine APRADD, 
the upper triangle of the matrix is used to form the a priori terms . In subroutine 
KALMAN , the SUM matrix is used for scratch storage . After the last iteration , the 
Cramer-Rao bounds and correlations are computed in the SUM matrix . 


1M0ECK :>UiC0»1 

/SMirnM/ 35, 35 ) 


SUMCOM* 1 

SUMCTK 2 
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A. 1.1. 46 


A . 1 . 1 . 46 SUMS AV . - Common deck SUMSAV contains information about the 
a priori penalty function . WAPR and ITAPR are described in reference 1 , sec- 
tion 3.3.8(22) . DIAGON is the weighting vector of squared elements selected from 
the APR prefix matrices (ref. 1, sec. 3.3.11(4)) . APRDIF is the vector of differences 
between the estimates and the a priori values. DIAGON is defined and APRDIF is 
initialized by subroutines ALLOW and VARY . 


CUftflOM /i'JHjAV/ oUGMi 35),4PrlOIF( ^5 ) . W»Pk. 


SUMSAV 1 

SUMSAV B 


A. 1.1. 47 TAPPOS . - Common deck TAPPOS contains information about the 
position of the time history data file. ITM is the last time read in total milliseconds; 
it is initialized to 0 in the main program . REW is used to request that subroutine 
READTH rewind the time history data file . REW is set to TRUE on the first point of 
a maneuver if the maneuver start time is less than or equal to the last time read . At 
all other times , REW will be FALSE . User subroutine READTH is responsible for 
checking REW and manipulating the data file as desired . 


♦CTMDECK TAPPOS TAPPOS t 

COMMON /TAPPOS/ ITM, PEW TAPPOS 2 

LOGTCAL PPW TAPPOS 3 


A. 1.1. 48 THPLOT$ . - Common deck THPLOT$ contains variables used only in 
subroutine THPLOT . If the program is run without segmentation or overlay , it 
may be desirable to shorten these vectors and store them in common blocks not used 
during the plotting (SUMCOM is the largest such block) . Time from the maneuver 
start is stored in the vector TIME . Measured and computed observations are stored 
in X and XX , respectively . NCH is the number of observation time histories stored 
simultaneously , and NTPLT is 2 plus the maximum number of time points plotted 
(see sec. 3.2) . XXX is equivalenced to X and XX; it is used to store up to 2 X NCH 
state, control, or extra signal time histories. Z , ZZ , and DC are used to read in 
each point of the time histories . VMINS and VMAXS contain the minimum and maxi- 
mum values of the signals plotted. IPLT is a vector used to indicate which of the 
states , controls , and extra signals are to be plotted . 


♦CTMDECK THPITTI 

OIMENSICN 2( 061,221 08),DC( 31),1PUT( 31I,VMIN»( 06I,VMAXS( 061, 
TlHE(1202),XXXn202, 06),X(1202, 03),XX(1202, 03) 

EOUIVALENCE ( X (1 , 1 ) > XXX ( 1 , 1 ) ) , ( XX (1, 1 ) , X XX ( 1, OA ) ) 

r. 

NTPLT«1202 
NCH" 03 


THPlOTi 1 
THPLOTS 2 
THPIOTS 3 
THPLOTS A 
THPLOTS 5 
THPLQTS 6 
THPLOTt T 


A . 1 . 1 . 49 TQDATA . - Common deck TODATA contains information needed to 
read the time history data file (channel numbers are in block INORD) . STC and ETC 
are the requested maneuver start and end times in total milliseconds . All of the other 
variables are input variables described in reference 1, section 3. 3. 8(3) , (8) , (9) , 
(29) , and (32) . The relative positions of ZBIAS , UBIAS , and EXBIAS in the common 
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block should not be changed , since subroutine THDATA depends on this relationship . 
The same applies to ZSCALE , USCALE , and EXSCAL. 


♦CnUDECK TODATA 

COMMON /TJOAFA/ STC ( 1 5 ) , E TC ( 15 ) * THIN# PR INTI> HAXRFC# 

ZBIAS( 08I#U9IAS( 0A)#EXRIAS( 20)#ZSCALE( 08)#U$CAL€( OA), 
exSCAK 20) 

INTEGER THIN#STC#ETC 
LOGItAl PKINTT 


TOOATA 1 
TODATA 2 
TOOATA 3 
TOOATA A 
TODATA 5 
TOOATA 6 


A. 1.1.50 TOGIRL . - Common deck TOGIRL contains input variables used to con- 
trol convergence . SNOISE is TRUE if the state noise algorithm is used; it is defined 
by subroutine ALLOW . The remaining variables are described in reference 1 , 
section 3.3.8(16) to (19) , and (21) . Variables that control G determination and 
a priori are in common blocks GICOM and SUMSAV , respectively . 


♦roMDECK TOGIRL TOGIRL 1 

COMMON /TOGIRL/ BOUND# ERRMAX# FUL LI #N0ITF R # OFAC # ITOFAC # SNOISE TOGIRL 2 

LOGICAL FULLl#5NrnSE TOGIRL 3 


A . 1.1. 51 TOGRAD . - Common deck TOGRAD contains system vectors used in 
computing the time histories and gradients . The suffix 1 indicates a value at the 
beginning of the sample interval; the suffix 2 indicates a value at the end of the 
sample interval; the suffix 12 indicates the average of the values at the beginning 
and end of the sample interval. Names without suffixes indicate the end of the 
sample interval. XT is the predicted state, XH the corrected state, U the control, V 
and W the known forcing functions in the state and observation equations , Y the 
predicted observation, Z the measured observation, ZMY the residual, and ZMYFLT 
the filtered residual . 


♦ C'^MOECK TOGRAD TOGRAD 1 

COMMON /TOGRAD/ XTK 07)#XT2( 07)#XH2( 07># TOGRAD 2 

\fH 07)*V2( 07),VI2( 07)#U1( 0A),U2( 0A).U12( 0A)f TOGRAD 3 

Y( 08)#2MY2< 08)#ZMYFLT( 08)#W( 08) TOGRAD A 


A. 1.1. 52 TOPLOT . - Common deck TOPLOT contains input variables used to 
control plot scales and signals plotted . All of the variables except for RATIO are 
described in reference 1, section 3.3.8(35), (36), (37), and (39) to (44). RATIO is 
PLTFAC/2 if INCH (ref. 1, sec. 3.3.8(38)) is TRUE and PLTFAC/2.54 if INCH is 
FALSE . The relative positions of XMAX , UMAX , and EXMAX should not be changed , 
since subroutine THPLOT depends on this relationship . The same is true for XMIN , 
UMIN, andEXMIN. 


♦C3MDECK TOPLOT 

CaMMuN /rOPLOT/ZMAX( 08)#ZMIN( 08)#XMAX( 07)»UMAX( 0A)*EXMAX( 20), 
XMIN( 07|#UMIN( OA), EXMINT 20)#XPL3T( 07 1 , NUPL Tt NEXPLTt 
T1MESC#RATI0»PAPER 
LOGICAL XPLOT 


TOPLOT 1 
TOPLOT 2 
TOPLOT 3 
TOPLOT A 
TOPLOT 5 
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A . 1 . 1 . 53 VARDEF$ . - Common deck VARDEF$ defines the values of the physical 
matrix dimensions. This common deck is called only in subroutine VARDEF. Matrix 
physical dimensions are discussed in section 3 and in reference 1 , section 2 . 


♦ nHOECK VARDEFS 
HKXt • 07 

MAXZ • 08 

^^4XU « 0^ 

M*X% ■ 04 

LEX - 20 

NT ■ 35 

HAXTV • 50 

MAXKV ■ 15 

HAXHRO • 36 
M4XSPT « 11 


VAROEFS 1 
VAROEFS 2 
VAROFFS 3 
VAROEFS 4 
VAROEFS 5 
VARDEFS 6 
VAROEFf 7 
VARDE^S 8 
VAR0EF4 9 
VARDEE$ 10 
VAROFFS 11 


A. 1 .1.54 XSUMS . - Common deck XSUMS contains the averages and standard 
deviations of the corrected states for the last iteration. Cards 37 to 39 and 116 
to 118 of subroutine GIRL accumulate the sums in XSUM and the sums of the squares 
in X2SUM . Then cards 36 to 42 of SUMOUT compute and print the averages in XSUM 
and standard deviations in X2SUM . 


♦ '•QHOECK XSUMS 

CJMMUN /XSUMS/ 07)fX2SUM( 07) 


XSUMS 1 

XSUMS 2 


A . 1 . 2 Standard Aircraft Routines 

The common decks of the standard aircraft routines are described below in 
alphabetical order . 

A. 1.2.1 FLCOND . - Common deck FLCOND contains variables describing the 
flight condition. All of the variables except G are input variables described in 
reference 1, section 4.3.3(10) and (15) to (21). G is the acceleration of gravity, 

o 2 

32 . 172 feet/ second or 9 .80665 meters/second , depending on METRIC (ref. 1, 
sec. 4. 3. 3(3)). 


♦C1MDECK FLCOND HCONP 1 

CONNON /FLCOND/ OB AR» V» THFTA» PHI » ALPHA. NACH, PAR*N,C6, G FLCOND 2 

REAL NACH FLCOND 3 


A . 1 . 2 . 2 GFDEFS . - Common deck GFDEFS contains the default values for F and 
GGI defined by subroutine ONCE . FLON and FLAT are the longitudinal and lateral- 
directional defaults for F . GGILON and GGILAT are the corresponding defaults 
for GGI . 


♦COMDfcCK GFDEFS GF0EF$ 1 

COMMCN /GFDEFS/ GG1LAT( 09 , 08)#G6ILGN< 09# 08)»FLAH 08# 07)# GFDEFS 2 

FL0N( 08» 07) GFDEFS 3 


A . 1 . 2 . 3 GRA V . - Common deck GRAY contains the derivatives of the gravity 
terms in the ct and p equations . DGDT is the derivative of the gravity term in a with 
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respect to 0 . DGDP is the derivative of the gravity term in p with respect to cp . 
These quantities are defined by subroutine MAKEL and subsequently used in both 
MAKEL and MAKEVW for the linearization of the gravity terms . 


♦CTMOFCK 6B»V 

COM, ITS /GRAV/ DGuT*DGDP 


GPAV 1 

GRAV 2 


A. 1.2. 4 INERTS. - Common deck INERTS contains aircraft mass and geometry 
data . All of the variables except MASS and WTCG are input variables described in 
reference 1, section 4. 3. 3 (2) , (4) , (5) , (6) , (8) , and (9) . MASS is the weight 
(ref. 1, sec. 4. 3. 3 (7)) divided by the acceleration of gravity. WTCG is the reference 
center of gravity of the predicted derivatives . If a predicted-derivative file is 
not used (ref. 1, sec. 3.3.4), WTCG is undefined. 


♦CTMDFCK INFRTS 

COMMOS /INERTS/ IX,IY,IZ, I X Z » I XE» M ASS . AR FA * CHORD. S P AN # WTCG . S HIFT 
REAL tX.lY.IZ.IXZ.lXEf MASS 
LOGICAL SHIFT 


INERTS 1 
INEPTS 2 
INERTS 3 
INEPTS 4 


A. 1.2. 5 INSTR . - Common deck INSTR contains instrument positions and cor- 
rections . All of the variables except DCGFT are input variables described in 
reference 1, section 4.3.3(11) to (14). If SHIFT is TRUE, DCGFT is the distance of 
the flight center of gravity forward of the reference center of gravity in feet or meters . 
If SHIFT is FALSE or if there are no predicted data, DCGFT is 0. 


♦rOMDECK INSTP 

COMMON /iNiTP/ KALF.KB.XALF,XB»XAN,XAX,XAY,YALF.YB.YAN#YAX#YAY, 
ZALF.ZP.ZAN.ZAX.Z AY. DCGFT 
REAL KALF.KB 


INSTR 1 
INSTP 2 
INSTP 3 
INSTP A 


A. 1.2, 6 LONLAT . - Common deck LONLAT contains the logical variables LONG 
and LATR, described in reference 1, section 4. 3. 3(1) . The program forces LATR to 
be .NOT. LONG. 


♦COMOFCK LONLAT 

COMMON /LONLAT/ LONG.LATP 
LOGICAL LONG.LATR 


LONLAT 1 

LONLAT 2 

LONLAT 3 


A. 1.2. 7 UVCOM . - Common deck UVCOM contains the vector UVAR described 
in reference 1 , section 4.3.3(22) . 


♦COMDECK UVCOM 

COMMON /UVCDM/ UVAP( 04) 
INTEGER UVAR 


UVCCM 1 

UVCOM ? 

UVCOM 3 
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A. 2 


The subroutines will be described in the order that they appear on the listing 
of supplement 1 . The main program and the subroutines of the basic program are 
described first . Then follow the general utility routines for matrix manipulation , 
plotting, and time conversion. The standard aircraft user routines are described 
next , followed by the EISPACK routines . The flow of most of the routines is so 
simple that flow charts would be superfluous . Flow charts are given for those few 
subroutines for which they are useful . Refer to section A . 1 for descriptions of the 
variables in common blocks . 


A . 2 . 1 Basic Program 

The major routines of the basic program are listed in the order of their use . 

Major routines are defined as the main program , including all subroutines called 
directly from the main program . Subroutine GIRL is so important that it is con- 
sidered to be a major routine even though it is called from subroutine NEWTON , 
instead of from the main program. Each major routine is followed by the associated 
minor routines , listed in alphabetical order . The subroutine descriptions all refer 
to the cards at which the subroutines are called. 

A. 2. 1.1 MMLE3 . - MMLE3 is the main program. It contains the only call to 
common deck HISTORY , which describes the modification history of the program . 

The program card defines the files and buffer sizes . The program card and the 
variable definitions in common deck VARDEF$ must be changed in order to change 
file numbers . Some systems may not allow a program card , in which case it 
should be deleted . 

The program first calls VARDEF to initialize file numbers and matrix physical 
dimensions . The variable ITM (last time read on the input time history file) is then 
initialized to 0 . The variable PLTOPN indicates whether the plot file has been 
opened; it is initialized to FALSE . 

Cards 37 to 43 read the syntax check card (ref. 1, sec. 3.3.1) and the user 
routines control card (ref. 1, sec. 3.3.2) and define the variable UMOD. The 
rest of the section is skipped if UMOD is FALSE . If UMOD is TRUE , subroutine ONCE 
is called for any user initialization , and the predicted-derivative control card 
(ref. 1, sec. 3.3.4) is read to define WTFILE. Depending on the predicted-derivative 
control card, user routine WTIN is then called to read predicted- derivative data. 

Cards 56 to 70 loop until all cases have been analyzed (indicated by the variable 
LAST returned from subroutine MTLOAD) . If the plot file was opened, subroutine 
PLOT is called after termination of the loop in order to close it . 

A . 2 . 1 . 2 VARDEF . - Subroutine VARDEF defines the variables describing the I/O 
file numbers and the matrix physical dimensions. Common deck VARDEF$ is included 
to define the basic matrix physical dimensions . Other dimensions are then computed 
from the basic ones . VARDEF calls subroutine GETSET to initialize the matrix 
routines and define the maximum physical dimension allowed . VARDEF is called 
at card MMLE3 . 31 . 
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Initialize . 
Define dimensions 
and file numbers 


Read syntax check card 
and user routines control card 


No User 

routines? . 


Do user initialization 


Read predicted' derivative 
control card 



Read input for case 


Yes Syntax 

check 

^Vi^option?^ 



Plot time histories 
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A. 2. 1.3 


A. 2. 1.3 TITLES . - Subroutine TITLES reads the title card for a case (ref. 1, 
sec. 3.3.6) . It also calls the DATE and TIME routines (sec. 2.6) to find the date 
and time of the run for identifying the printout . Subroutine HEAD is then called to 
print the page heading. TITLES is called at card MMLE3.56. 

A . 2 . 1 . 3 . 1 HEAD . - Subroutine HEAD prints a page heading , consisting of the 
title , date , and time . HEAD is called from several different routines . 

A . 2 . 1 .4 EDIT (WTFILE) . - Subroutine EDIT reads the NAMELIST INPUT , the 
signal labels, and the time cards (ref. 1, secs. 3.3.8 to 3.3.10) . The first section 
of code defines basic program default values . The user routine USERIN is called at 
card 117 to read any input for the user routines (ref. 1, sec. 3.3.7) . The argument, 
WTFILE , is passed to subroutine USERIN to inform USERIN whether predicted- 
derivative data are available. USERIN may modify the basic program defaults. The 
next section of code reads the NAMELIST INPUT and makes some consistency checks 
between the options . The next sections print out scalar variables and options , read 
signal labels , and print vector variables and options . The last section of code reads 
and prints the requested maneuver times. EDIT is called at card MMLE3.57. 

A. 2 . 1 .5 MATSET . - Subroutine MATSET initializes the input matrices as re- 
quired by the standard matrix routines and defines their defaults . The N and V 
suffix and APR prefix matrices are initialized to 0 (except for RN , which is initialized 
to identity) . Each element of the M suffix matrices is initialized to 1. The GGI 
matrix is initialized to 0, and the hard constraint (HARD) and soft constraint (SOFT) 
matrices are initialized to indicate no constraints . The matrix labels and input 
flags in common block MATLAB are also defined by MATSET . MATSET is called 
at card MMLE3.58. 

A . 2 . 1 . 5 . 1 MTSET (AN , AV , APRA , AM , IM , II , JJ) . - Subroutine MTSET is 
used by MATSET to initialize a group of related matrices . IM is a code indicating 
which group of matrices is being initialized . II and JJ are the logical matrix dimen- 
sions to be used. For the matrices related to F (IM equals 9) , an FM matrix is not 
defined , so the code initializing the M suffix matrices is skipped . 

A . 2 . 1 . 6 MTLOAD (LAST) . - Subroutine MTLOAD controls the matrix input for 
each case . It reads the matrix header cards and determines which matrix is being 
read . It then sets the input flag (INFLAG) for that matrix and calls MATLD or 
CONIN to read the matrix body into the appropriate locations . Common block INMAT 
is used to pass information from the header card to MATLD and CONIN . Subroutine 
MTLOAD also detects the endcase card (ref. 1, sec. 3.3.12) , which signals the 
end of the matrix input . The variable LAST is defined based on the endcase card 
and passed back to the main program as an argument . This variable flags the last 
case of a run. MTLOAD is called at card MMLE3.59. 

A. 2. 1.6.1 CONIN (CON) . - Subroutine CONIN reads the body of a constraint 
matrix input into the argument , CON , and prints out the matrix . Depending on the 
matrix header card , the matrix input flag (INFLAG) is reset to 0 , allowing the 
constraints read in to supplement rather than replace any default constraints . This 
is discussed in section A. 2. 3. 6 and in reference 1, section 3.3.11. CONIN is 
called at cards MTLOAD . 85 and 87 . 
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A. 2. 1.6. 2 


A . 2 . 1 . 6 . 2 LOADED (AN) . - Logical function LOADED determines whether a 
given matrix has been read in from cards . The function LOADED is intended for 
use by user routine MATDEF to determine whether the matrix defaults are used . The 
matrix itself is used as an argument , and TRUE is returned if the matrix was read 
from cards. LOADED first calls GETLAB to extract the matrix name from the matrix, 
and calls MATNO to find the corresponding matrix number . The vector of input flags 
(INFLAG) is then checked to see if that matrix was read . LOADED will return the 
value FALSE for constraint matrices that should supplement any default constraints 
rather than replace them . This is because subroutine CONIN has reset the corre- 
sponding input flag to 0 . LOADED is called many times in the standard aircraft 
routines MATDEF and WTDEF . 

A. 2. 1.6. 3 MATNO (ALAB) . - Function MATNO returns a matrix number, given 
its name as an argument . MATNO searches the list of names in common block 
MATLAB to find a matching name . The value returned is the index of the matching 
name . An error message is printed if the name is not found in the list . MATNO is 
called from several different routines . 

A . 2 . 1 . 7 THDATA . - Subroutine THDATA controls reading and processing of 
the input time history file UDATA . Subroutine READTH is called to do the actual 
manipulation of the data file so that the input format can be easily changed . THDATA 
handles the time searching, data scaling, printing, averaging, and associated tasks. 
For each maneuver (ref. 1, sec. 3. 3. 8(1)), THDATA defines the variable REW based 
on the last time point read and the requested start time . It then enters a search 
loop from cards 48 to 55 for the start time. The actual start time used for each 
maneuver is printed and stored for use by subroutine THPLOT . 

Cards 58 to 91 then process data until the maneuver stop time is found. The 
data are thinned if desired . At each time point , the requested data channels are 
extracted from the input record , and scale factors and biases are applied . Subrou- 
tine THMOD is called to modify or correct the data as desired . The data are then 
written to the scratch file UTl and printed if desired. Based on the first and 
second thinned time points of the first maneuver , THDATA computes and prints the 
sampling rate of the data as described in reference 1, section 3. 3. 8 (4) . 

After all of the maneuvers have been read, THDATA computes and prints the 
averages over the entire case of the observations, controls, and extra signals. 
Standard deviations, minima, and maxima are also computed, but not printed. User 
routine AVERAG is then called to allow the user access to these averages. 

Subroutine THDATA stops with an error message if no time points are found in 
a requested interval . Another error check limits the number of calls to user routine 
READTH for each case; this is to guard against possible infinite loops caused by 
logic errors or omissions in READTH . End-of-file checks are made if user routine 
READTH has defined the variable EOFTH . End-of-file is not considered an error in 
itself, but may result in an error if no time points are found in an interval or if a 
following interval is requested . 

THDATA is called at card MMLE3.61. 

A . 2 . 1 . 8 COMP AT . - Subroutine COMP AT sets the logical matrix sizes to com- 
patible values and checks dimension limits. Cards 24 to 48 determine the matrix sizes 
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to be used as described in reference 1, section 3.3.8(11) to (14) and check these 
sizes against the dimension limits. Cards 50 to 78 set the appropriate logical sizes 
to be used for each of the matrices . Card 79 calls SYM to check GGI for symmetry . 
The starting nondimensional matrices and the F and GGI matrices are then printed 
out. COMP AT is called at card MMLE3.63. 
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A ■ 2 . 1 . 9 ALLOW . - Subroutine ALLOW determines what parameters are allowed 
to vary. It also defines maps between different subsets of the parameters . This 
information is obtained from the V suffix matrices and constraint matrices . It is 
output from ALLOW in the DETERM and MAPCOM common blocks . ALLOW can be 
viewed as translating information from the input format (organized for ease of use) 
to the internal format (organized for compactness and efficiency) . 

ALLOW first initializes several vectors. It then calls VARY to define the 
independent unknowns in the N suffix matrices and the F matrix. Cards 41 to 50 
add variable initial conditions to the list of independent unknowns as requested by 
the input vector VARIC . ALLOW then checks the dimension limit on the number of 
independent unknowns . 

Subroutine CONSTR is called for preliminary processing of the hard and soft 
constraint matrices . Hard constraint processing is then completed by the call to 
HARDC . The remainder of the soft constraint processing is done later in subroutine 
APRADD . 

Cards 59 to 63 determine whether the state noise algorithm will be used. Then 
cards 71 to 91 define the maps in common block MAPCOM that are used in the state 
noise algorithm. Cards 93 to 97 print the results from subroutine ALLOW. 

ALLOW is called at card MMLE3.64. 

A. 2. 1.9.1 CONSTR (CON) . - Subroutine CONSTR does preliminary processing 
of the soft or hard constraint matrices . The matrix to be processed is the argument . 
The matrix names read in to define the constraints are translated into matrix numbers 
(see sec . A . 1 . 1 . 8) by calls to subroutine MATNO . If constraint ratios were not 
specified on input , they are defined from the ratios of the starting values . Error 
messages are provided for unallowed matrix numbers and ill-defined constraint 
ratios. CONSTR is called at cards ALLOW. 55 and 57. 

A. 2. 1.9. 2 GVALVE (IM , IR , IC) . - Function GVALVE returns the present value 
of any coefficient in the N suffix or F matrices . Input arguments are the matrix num- 
ber (see sec. A. 1.1. 8), row, and column (IM , IR, and IC, respectively). GVALVE 
is called at cards CONSTR. 27 and 32. 

A . 2 . 1 . 9 . 3 HARDC . - Subroutine HARDC implements the hard constraints . 

Before HARDC is called , the lists in common block DETERM must have been defined 
for the independent unknowns . Subroutine CONSTR must have been called for 
preliminary processing of the hard constraint matrix. HARDC extends the list in 
common block DETERM to include the hard constraints . For each hard constraint , 
subroutine LOCATE is called to locate the independent variable of the constraint in 
common block DETERM. If the independent variable is found, the constraint infor- 
mation is added to the lists in DETERM . The variable NVAR , specifying the length 
of the lists in DETERM , is set to the number of independent unknowns during the 
execution of HARDC . This is so that the search in subroutine LOCATE will be 
restricted to the independent unknowns . At the end of subroutine HARDC , NVAR is 
set to the number of independent unknowns plus the number of active hard con- 
straints . An error message is provided for exceeding the dimension limits in 
common block DETERM . HARDC is called at card ALLOW .56. 
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A. 2. 1.9. 4 LOCATE (IM, IR, IC) . - Function LOCATE returns the index of a 
variable in the list of unknowns . The input arguments are the matrix number , row , 
and column (IM , IR , and IC , respectively) . If the given variable is not found in 
the lists in common block DETERM , the value 0 is returned . LOCATE is called at 
cards HARDC . 19 and APRADD . 20 and 21 . 

A. 2. 1.9. 5 SETCON (CON, AI , IR, IC , AJ , JR , JC , FACT, IFZERO) . - Subroutine 
SETCON defines a single default constraint in the hard or soft constraint matrices . 
SETCON i's intended for use in user routine MATDEF . The first argument is the 
hard or soft constraint matrix (HARD or SOFT) . AI , IR , and IC specify the location 
of the dependent variable in terms of matrix, row , and column, respectively. For 
convenience of use , the matrix itself (rather than the matrix name or number) is 
used for the second argument. Similarly , AJ , JR, and JC specify the matrix, row , 
and column of the independent variable . FACT is the constraint ratio . The last 
argument , IFZERO , is relevant only if the constraint ratio given is 0 . If the con- 
straint ratio given is 0 and IFZERO is FALSE , the constraint will be ignored as 
irrelevant . This is the outcome usually desired if the constraint ratio is a calculated 
quantity which can validly be 0 . If the constraint ratio given is 0 and IFZERO is 
TRUE , the constraint will be retained . Subsequent processing by subroutine 
CONSTR will define the constraint ratio as the ratio of the starting values . SETCON 
is called several times in the standard aircraft routine MATDEF . 

A. 2. 1.9. 6 VARY (AN, AV, APRA, AR) . - Subroutine VARY determines what 
coefficients in a given matrix are independently varying . The four arguments are 
the nondimensional starting matrix (AN) and its associated variation (AV) , a priori 
weighting (APRA) , and a priori value (AR) matrices . The corresponding V suffix 
matrix is searched for nonzero elements . For each such element , the matrix number , 
row, and column are stored in common block DETERM; corresponding a priori 
information is stored in common block SUMSAV . If there are no independent un- 
knowns in the matrix, the matrix name will be defined as DONT; this will prevent 
printing the matrix every iteration. Subroutine VARY distinguishes between the 
starting values (AN) and the a priori values (AR) , although the MMLE3 program 
does not currently preserve this distinction. VARY is called at cards ALLOW. 32 to 40. 

A . 2 . 1 . 10 NEWTON . - Subroutine NEWTON controls the iteration for obtaining 
the maximum likelihood estimates. NEWTON is called at card MMLE3.65. The 
following are important variables used in the iteration control: 

Variable Name Description 

ITA Number of iterations remaining until a priori is turned off plus one . 

When ITA is 0, a priori will either remain on or is already off. 

DFAC multiplication (ref. 1, sec. 3.3.8(21)) and G determination 
(ref. 1, sec. 3.3.8(24)) cannot start until ITA is 0. If conver- 
gence is achieved while ITA is nonzero , then ITA is set to 0 , 
a priori is turned off, and iteration continues. 

ITD Number of iterations remaining with DFAC multiplication . DFAC 

multiplication does not start until ITA is 0 . DFAC multiplication 
is not done the first iteration , regardless of ITA , unless FULLl 
is TRUE. The convergence test is disabled while DFAC is used. 

G determination cannot start until ITD is 0 . 
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IXGGI Number of iterations until G determination starts. If ITGGI is 0, 

either G determination has already started or will not be used 
(depending on ITG) . ITGGI does not start counting iterations 
until ITA and ITD are 0 . If convergence is achieved while ITGGI 
is nonzero and ITA is 0 , then ITGGI is set to 0 , and G determina- 
tion starts. 

CONVRG Convergence indication . CONVRG is set to TRUE when the cost 

functional converges within the limit specified by BOUND . 

CONVRG is used to turn off a priori , turn on G determination , or 
stop iteration , depending on ITA and ITGGI . 

The iteration loop is skipped by cards 23 to 29 if NOITER is 0 or if there are no 
unknowns. Cards 30 to 33 initialize iteration control variables. Cards 35 to 92 
are the iteration loop . 

GIRL is called at card 37 to compute the time history and the gradients of the 
cost functional. RESIDE is then called at card 39 to compute and print the residual 
powers and related quantities . Cards 42 to 44 determine if convergence has occurred; 
this determination may subsequently be changed, depending on the options in effect. 

Cards 46 to 53 control the a priori option. If ITA (initialized at cards 30 and 31) 
is nonzero , it is decremented by 1 each iteration . When ITA reaches 0 or convergence 
occurs, a priori is turned off by setting WAPR and ITA to 0; the convergence flag 
is turned off so that iteration can proceed with a priori off. On subsequent iterations, 
ITA is 0, so this logic is skipped and a priori remains off. If ITA is initially 0 , the 
code to turn off a priori will never be executed; the a priori weighting will thus 
remain at its initial value , which could be zero or nonzero . 

Cards 55 and 56 control the call to BIAS to determine linear unknowns only. 

Note that the definitions of the convergence flag CONVRG on cards 42 and 86 check 
FULLl or FULLIT to insure that an iteration on which BIAS is called will never be 
judged to have converged. Also card 58 checks FULLIT so that the call to DFACT 
will be skipped on iterations that BIAS is called . Finally , the call to FLIMIT at 
card 68 is not needed on iterations with BIAS because none of the linear unknowns 
affect the Kalman gain . 

Cards 58 to 62 control the call to DFACT, which implements the diagonal con- 
vergence factor option . DFACT is not called until ITA is 0 . DFACT is also not 
called on iterations for which BIAS was called . For each iteration that DFACT is 
called, ITD is decremented by 1; the convergence flag is forced to FALSE because 
convergence with DFACT is very slow and would often set the convergence flag 
before truly converging. After ITD reaches 0 (or if it starts at 0) , subsequent 
iterations do not call DFACT . 

Cards 64 and 65 invert the second gradient matrix and fill in its symmetric, 
upper triangular portion. Then cards 66 and 67 compute and print the changes in 
the parameter estimates . Card 68 calls FLIMIT if required to implement the inequality 
constraints described in reference 1, section 1.2.3. Cards 69 and 70 then revise 
and print the parameter estimates . 
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Cards 72 to 86 do the G determination if requested by a nonzero value of ITG. 

This code is not entered until both ITA and ITD are 0 . After ITA and ITD are 0 , 
cards 74 to 76 decrement ITGGI by 1 each iteration until ITGGI reaches 0 or converg- 
ence is obtained; either of these conditions triggers the start of G determination . 

The call to GIRL at card 82 computes the time history (no gradients are computed at 
this call) . RESIDS is then called at card 85 to revise GGI based on the residuals. 

The convergence is tested at card 86 . When G determination is active , each iteration 
has two steps: First, cards 37 to 70 revise all of the estimates except for GGI; and 
second, cards 82 and 85 revise the estimate of GGI. 

Cards 89 to 91 exit the iteration loop if final convergence has been attained. 

Card 94 prints a warning message if the iteration limit is reached without attaining 
convergence . 

GIRL is called at card 98 to compute the final time histoiry and gradients (the 
second gradient will be required to compute the Cramer-Rao bounds) . RESIDS , 
called at card 100 , computes and prints the final iteration residual powers and re- 
lated quantities. 

A . 2 . 1 . 10 . 1 APRADD . - Subroutine APRADD adds a priori terms to the first and 
second gradients . Soft constraints are implemented by APRADD as off-diagonal 
a priori terms . The upper triangular part of the SUM matrix is used to form the 
a priori weighting matrix; the diagonal elements of this matrix are stored in row JKM 
of SUM . APRADD is called at card NEWTON .53 . 

A . 2 . 1 , 10 . 2 BIAS . - Subroutine BIAS causes only bias and control terms (linear 
terms) to be estimated in a particular iteration. The logical variable FULLIT controls 
the call to BIAS at card NEWTON .56 . FULLIT is defined in turn at card NEWTON .55 
depending on FULLl and the iteration number. This results in BIAS being called 
for the first iteration unless FULLl is TRUE . 


A . 2 . 1 . 10 . 3 DFACT . - Subroutine DFACT implements the diagonal convergence 
factor option (ref. 1, sec. 3.3.8(21)) . It multiplies the diagonal elements of the 
second gradient by DFAC . The calling of subroutine DFACT is controlled by the 
input variable ITDFAC (ref. 1, sec. 3.3.8(21)). Cards NEWTON .58 to 62 implement 
the logic to call DFACT . 


A . 2 . 1 , 10 . 4 FADJ (AVG) . - Subroutine FADJ adjusts F during G determination . 
The intent is to keep the Kalman gain matrix, K, unchanged as closely as reasonable. 
On entry , it is assumed that WFRSQ contains the diagonal elements of .^GGF^^^TggF^^ 

and that the diagonal FRSQ contains the old GGI matrix elements; both of these quan- 
tities were computed in subroutine RESIDS. The input argument, AVG, is assumed 
to contain the logarithmic average of the elements in WFRSQ . The algorithm used is 
to multiply each row of F by the ratio of the corresponding diagonal element of 


ERIAC* GGI^j^ DIAG (WFRSQ) ERIAC 

ERIAC* GGI , , ERIAC 
old 


where only diagonal elements of GGI^^^ 


are used 


and DIAG (WFRSQ) is the diagonal matrix formed from WFRSQ. If any element of the 
above numerator is 0 , the corresponding row of F is instead multiplied by AVG. Only 
independently varying elements of F will be changed by subroutine FADJ as controlled 
by the loop from cards 33 to 41. FADJ is called at card RESIDS .95. 
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A. 2. 1.10. 5 FLIMIT . - Subroutine FLIMIT constrains certain diagonal elements 

of the closed loop gain K(ER + C) to be less than or equal to 1. The algorithm 
and reasons for this constraint are discussed in reference 1, sections 1.1.2 and 1.2.3. 
GRADX, GRADY, and GRADl are used for scratch storage in this routine. The 
subroutine first forms the gradients of these diagonal elements in GRADX. It then 
determines which diagonal elements of F are varying, using DUM2X to flag such 

elements . Constraints are only made on diagonal elements of K {ER ^A + C) corre- 
ponding to unknown diagonal elements of F. Cards 64 to 76 compute the linearized 

extrapolation of the diagonal elements of K(,ER~^A + C) , adding the current value 
plus the gradient times PB (the vector of proposed coefficient changes) . The rows 
of GRADX corresponding to constraints that are satisfied (or the elements that are not 
constrained) are deleted, and the remaining rows are compressed into the first JJ 
rows of GRADX. Corresponding elements of DUMX are filled with the amount by which 
the constraint is exceeded. If no constraints are exceeded (JJ equals 0) the sub- 
routine is done. Otherwise, cards 79 to 89 modify PB to lie on the constraint boundary , 
approximated by local linearization. FLIMIT is called at card NEWTON. 68. 

A. 2. 1.10. 6 RESIDE (GUT) . - Subroutine RESIDE does computations based on the 
sample residual power . It first computes the filtered and unfiltered sample residual 
powers from the accumulated sums at cards 21 to 27. Cards 28 to 37 then eliminate 
the effect of any unweighted signals; this is needed so that such residuals do not 
affect the inverse and determinant of the residual power. Next, cards 40 to 54 
compute and print the weighted errors for the filtered and unfiltered residuals . The 
cost functional (unfiltered weighted error sum) is placed in ERRVEC (NITER) and 
the filtered weighted error sum is placed in ERRFLT . WRSQ and WFRSQ are the 
diagonal elements of the unfiltered and filtered weighted errors , respectively . 

Cards 56 to 62 compute and print the log determinant of the unfiltered residual power 
using the product of the eigenvalues . 

The remainder of the subroutine does the G determination . The subroutine 
argument , GIIT , determines at card 66 whether this code is bypassed . The old GGI 
matrix is saved in FRSQ. Cards 68 to 72 move the unfiltered sample residual power 
matrix (or its diagonal elements only) into the GGI matrix and then invert it. This 
is the preliminary value of the new GGI matrix. Cards 73 to 92 then apply relaxation 
to further revise GGI; cards 80 to 91 compute WFRSQ as will be required by subroutine 
FADJ . After symmetrizing and printing the new GGI matrix at cards 93 and 94 , sub- 
routine FADJ is called at card 95 if the state noise algorithm is used . FADJ will 
adjust F to compensate for the GGI change. 

RESIDE is called at cards NEWTON . 39 , 85 , and 100 . 

A. 2. 1.10. 7 SPITEM. - Subroutine SPITEM prints the N suffix matrices and the F 
matrix . Subroutine VARY controls which matrices are printed by changing some of 
the matrix names to DONT . Subroutine SPIT will ignore any matrices which have 
DONT as a name. The call to SPITEM at COMPAT.82 occurs before VARY is called; 
therefore , all of the matrices will still have their proper names and thus will be 
printed. The call to SPITEM at NEWTON. 70 is subsequent to VARY; therefore, some 
(or all) of the printouts may be omitted . 
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A . 2 . 1 . 10 . 8 UPDATE . - Subroutine UPDATE updates the parameter estimates . PB 
contains the vector of parameter changes to be made. Cards 43 and 44 keep track of 
the total change from the a priori value for use in subroutine APRADD . UPDATE 
is called at card NEWTON. 69. 

A. 2. 1.11 GIRL (DOGRAD, LAS TIT ) . - Subroutine GIRL computes time histories 
and , optionally , gradients . The argument DOGRAD controls whether gradients are 
computed; the argument LASTIT controls various output options used only on the 
last iteration. Cards 31 to 63 do initialization before entering the case and time loops. 
Average dimensional matrices are computed by REAT or DIM2 if needed. If the time 
varying option is not used or if test output is requested , the average dimensional 
matrices and transition matrices are required, so REAT is called. If REAT is not 
called and the state noise algorithm is used, DIM 2 is called to compute the matrices 
required for computation of the Kalman gain matrix and its gradient by subroutines 
KALMAN and GRADK . Cards 50 to 56 test if GGI is diagonal . Cards 66 to 76 initialize 
each maneuver and do any output required for the first time point of the maneuver . 
Subroutine INIT is called to define the time history initial condition , and GRADIC 
defines the gradient initial condition . 

The time loop goes from cards 78 to 134. The measured data are read from 
scratch file UTl and, if TIMVAR is TRUE , REAT is called to recompute the dimen- 
sional matrices at each time point . The predicted response is computed at cards 82 
to 95 . The filtered and unfiltered residuals are then computed and summed . The 
unfiltered residual is also stored as an augmented column to the GRADY matrix for 
convenience in computing the cost functional gradients . If the state noise algorithm 
is used , the corrected responses are then computed . The gradients are computed at 
cards 123 to 129. Subroutine GRAD computes the gradient of the predicted response 
at one time point and stores it in the matrix GRADY . The call to SUMULT at card 129 
accumulates the contribution from the time point to the first and second gradients of 
the cost functional . 

After the case and time loops , cards 137 and 138 move the first gradient of the 
cost functional from row JKM to column JKM of the SUM matrix , as required by the 
rest of the program . DIM2 is called to recompute the average dimensional matrices 
if required for subroutine FLIMIT . 

GIRL is called by cards NEWTON .37, 82, and 98. 

A. 2. 1.1 1.1 CALLAM. - Subroutine CALLAM stores the dimensionalization ratios 
in common block AMCOM . It calls user routine MAKEM if appropriate to compute 
the M suffix matrices. CALLAM is called from card DIM1.18. 

A. 2. 1.11. 2 DIMl. - Subroutine DIMl computes the basic dimensional system 
matrices . If USEAVG is TRUE , it moves average values into common block BILIN so 
that average dimensional matrices will be computed . CALLAM computes the M suffix 
matrices, and MAKEL adds in the L suffix matrix contributions. DIMl is called from 
card DIM2.11. 

A. 2. 1.1 1.3 DIM2. - Subroutine DIM2 computes dimensional system matrices and 
expressions . It calls DIMl to compute the dimensional matrices . All of the matrix 
expressions in the common blocks ECOM and OBSRV are then evaluated . DIM2 is 
called from cards REAT. 13, INIT. 26, and GIRL. 47. 
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A . 2 . 1 . 11 . 4 GRAD . - Subroutine GRAD computes the gradient of the predicted 
response for one time point. Cards 23 to 26 compute the variables in common deck 
GRAD$, used later in the subroutine. The gradients of the predicted response are 
then computed using the equations given in reference 1, section 3.1. It is assumed 
that, on entry , GRADX contains the gradient of the corrected response at the previous 
time point. Cards 30 to 34 zero GRADY and multiply GRADX by the transition matrix. 
GRADl is used for scratch storage. The loop from cards 35 to 66 adds terms to 

GRADX and GRADY for each unknown . Then card 67 adds the (ER + C) GRADX 
term into GRADY . 

At this point GRADX and GRADY contain the gradients of the predicted state and 
observation. If the state noise algorithm is used, cards 72 to 81 add 
-K Vz + (VK) (z - z) to the gradient of the predicted x to obtain the gradient of the 
corrected x. This is needed for the next call to subroutine GRAD; only GRADY is 
required as an external output from GRAD . If the state noise algorithm is not used, 
the corrected and predicted states are identical , so the last section of code is skipped . 

GRAD is called from card GIRL. 123. 

A. 2. 1.11. 5 GRADIC . - Subroutine GRADIC initializes the gradient of x at the 
beginning of each maneuver . All elements of the gradient are 0 , except those 
corresponding to variable initial conditions. GRADIC is called from card GIRL. 73. 

A . 2 . 1 . 11 . 6 GRADK . - Subroutine GRADK computes the gradient of the Kalman 
gain matrix. The results are stored in the triply dimensioned array DK in common 
block GRDCOM . GRADK assumes that subroutine KALMAN was previously called to 
compute the Riccati covariance matrix, P, and the Kalman gain matrix, KGAIN. 

GRADK first computes RIAP as RIA X P , and calls GRADP . RIAP is used both in 
GRADP and GRADK . On return from GRADP , the gradient of P is stored in DK . 

Cards 25 to 28 multiply VP by CTG = (ER + C)*GGI (computed in subroutine 

KALMAN) . Cards 30 to 54 then add in the elements of VP(ER + C)*GGI to com- 
plete the gradient of K. Finally , DK is printed if the TEST option is on. 

GRADK is called from card GIRL. 57. 

A . 2 . 1 . 11 . 7 GRADP . - Subroutine GRADP computes the gradient of the Riccati 
covariance matrix and stores the result in the array DK . It assumes that the matrices 

RIAP = R~^AP, RIF = R~^F, and RIFRIF = R~^F(R~^F)* have been computed by the 
subroutines GRADK and KALMAN . RIFRIF will be destroyed in GRADP . As discussed 
in reference 1 , section 1.6, the gradient of the Riccati covariance matrix is computed 
as the solution to a group of Lyapunov equations . 

Cards 24 to 58 compute one-half of the constant terms in the Lyapunov equations 
and store the results in DK . To get the full constant terms , the values stored in DK 
would be added to their transposes . Cards 60 and 61 compute the coefficient of VP in 
the Lyapunov equations and store the result in DUM. Cards 65 to 74 perform the real 
eigenvector decomposition of DUM . RIFRIF becomes the block-diagonalized DUM 
matrix, DUM2 is the matrix of eigenvectors transposed, and DUM3 is the inverse 
of DUM2 . 
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The loop from cards 76 to 88 solves one of the Lyapunov equations for each pass 
through the loop . Cards 77 to 82 compute the constant term of the block-diagonalized 
Lyapunov equation . Then LYAPCB solves this block-diagonalized equation . Cards 84 
to 87 transform this solution back into the solution of the original Lyapunov equation . 

GRADP is called from card GRADK . 23 . 

A . 2 . 1 . 11 . 8 INIT (IT , XT 1 , U1 , Y , V , W) . - Subroutine INIT sets the time 
history iriitial conditions and biases. The arguments, IT, XTl, Ul, Y, V, and W, 
are the total time in milliseconds , state , biased control , computed observation , state 
equation forcing function , and observation equation forcing function , respectively , 
all at the initial time point of a maneuver . Card 16 reads the measured data and stores 
them in common block BILIN. Cards 18 to 25 define the default state, observation 
bias , and control appropriately for perturbation equations . If TIMVAR is TRUE , 

DIM2 is called to compute dimensional matrices at the initial time point . The V 
and W vectors are computed for use in defining the computed observations. Card 28 
calls subroutine UINIT to redefine the state , observation bias , and control as desired 
by the user. Cards 30 to 32 compute the control bias as the difference between the 
model control and the measured control at the initial point. Variable initial condition 
increments are added and OBSERV is called to compute the initial computed observation . 
INIT is called from card GIRL. 68. 

A. 2. 1.1 1.9 KALMAN. - Subroutine KALMAN computes the Kalman gain matrix. 

The matrices RIF , RIFRIF , and CTG computed here are used subsequently in sub- 
routines GRADK and GRADP . Cards 20 to 22 define the physical dimensions of P, 

KGAIN, and CTG. Then cards 23 to 28 compute the matrices of the continuous-time 
Riccati equation as follows: 


RIFRIF = R V(R ^F)* 

DUM = (ER~^A + C)*(GG*)~^(ER~^A + C) ^ 

Cards 32 to 35 destroy the physical dimension pointer of SUM and define two scratch 
matrices in the space occupied by the SUM matrix . RICATC is then called to solve 
the continuous Riccati equation. DUM2 and the two scratch matrices defined by SUM 
are used for scratch storage matrices in RICATC . The last two arguments of the 
call to RICATC are columns of SUM used as scratch vectors in RICATC . The scratch 
usage of the SUM matrix here assumes that the physical number of rows of the SUM 
matrix is at least two times the physical dimension of the state vector plus one, i.e. , 
NI = 2 X MAXX + 1 (see sec. 3) . Card 43 multiplies the Riccati covariance matrix 

by (£R + C)*(GG*) ^ to obtain the Kalman gain matrix. 

KALMAN is called at card GIRL . 48 . 

A . 2 . 1 . 11 . 10 OBSERV (X , U , ONES , V , W , Y) . - Subroutine OBSERV expands 
the observation equation. X, U, ONES, V, and W are input arguments giving the 
state , control , bias , state equation forcing function , and observation equation 
forcing function , respectively . The result is stored in Y . OBSERV is called at 
GIRL . 96 and INIT . 36 . 
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A . 2 . 1 . 11 . 11 RE AT . - Subroutine RE AT computes dimensional matrices and 
transition matrices . It first calls subroutine DIM2 to compute the basic dimensional 
system matrices and matrix expressions . Then EAT is called to compute the transition 
matrix, PHI, and its integral, BUM. PSI, PSIB, and PSIS are expressions involving 
the integral of the transition matrix . REAT is called at cards GIRL . 44 and 80 . 

A . 2 . 1 . 11 . 12 SPIDIM . - Subroutine SPIDIM prints the dimensional system matrices , 
matrix expressions , and transition matrices . Subroutine REAT must be called before 
SPIDIM in order to define all of the matrices printed . SPIDIM is called at card GIRL . 45 . 

A. 2. 1.12 SUMOUT . - Subroutine SUMOUT handles various output to summarize 
a case (not including plotted output) . The variable BLOWUP , defined by subroutines 
NEWTON and GIRL , is first examined to see if the case has diverged past program 
limits. If so, most summary output is skipped, plotting is turned off, and measured 
time histories are printed if desired . The rest of the subroutine is only executed if 
BLOWUP is FALSE . 

Subroutine CRAMER is called at card 18 to compute Cramer-Rao bounds and 
parameter correlations . OUTRUN is called at card 19 , depending on PUNCH , to 
punch coefficient estimates as required . Both CRAMER and OUTRUN are bypassed 
if NOITER is 0. Cards 21 to 34 normalize the residual covariance matrix to obtain 
the residual correlations and print them out . The unnormalized covariance matrix 
is saved in FRSQ during this computation . Cards 36 to 42 compute and print the 
averages and standard deviations of the corrected state estimates . The summary of 
the convergence of the cost functional is then printed . Finally , the plotting error 
limit is checked; if it is exceeded , plotting is turned off and , optionally , measured 
time histories are printed . 

SUMOUT is called from card MMLE3.66. 

A . 2 . 1 . 12 . 1 CRAMER . - Subroutine CRAMER computes the Cramer-Rao bounds 
and estimated parameter correlations . The Fisher information matrix must be in 
SUM when CRAMER is called . First the matrices in common block CRMAT are 
initialized to 0 . Then the SUM matrix is inverted. The loop from cards 28 to 53 
places the Cramer-Rao bounds in the appropriate locations in the matrices in CRMAT . 
This loop assumes that the independent unknowns are the first JKMMl unknowns in 
the list in common block DETERM. These matrices are then printed. Cards 64 to 75 
normalize the inverse of the information matrix to obtain the estimated correlations , 
which are then printed . CRAMER is called from card SUMOUT . 18 . 

A . 2 . 1 . 12 . 2 CRSET (AC , AN , ALAB) . - Subroutine CRSET initializes one of the 
matrices of Cramer-Rao bounds . The argument ALAB is used for the name of the 
matrix AC unless the name of AN is DONT . If the name of AN is DONT , the name of AC 
is also set to DONT; this will prevent AC from being printed. CRSET is called at 
cards CRAMER. 16 to 24. 

A . 2 . 1 . 12 . 3 ERRTHP . - Subroutine ERRTHP prints the measured time history . 

It is called at card SUMOUT .48 if the cost functional exceeds allowable limits. The 
data printed are from scratch file UTl, which has already had scale factors, biases, 
and modifications from subroutine THMOD applied . 
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A .2 .1.13 THPLOT (PLTOPN) . - Subroutine THPLOT plots time history data. 

The data to be plotted are obtained from the scratch file UT2 . The argument, 

PLTOPN , controls whether the plot file is to be opened by the call to PLOTS at 
card 23. If PLTOPN is TRUE , the file is already open, so the call to PLOTS is 
skipped . Cards 34 to 44 put the title , start time , maneuver number , and other iden- 
tifying information on the plot . Cards 46 to 56 determine the time scale and thinning 
used. Then cards 59 to 65 form the time vector and draw the time axis and label. 

The 'observations are read into core and plotted up to NCH at a time by cards 67 
to 110. Cards 68 and 69 determine how many observations are left to plot (if any) . 
Then cards 71 to 74 position the file UT2 at the beginning of the maneuver being 
plotted . (NIP is the total number of time points in the previous maneuvers . ) 

Cards 76 to 88 read NCHAN measured and computed observations into the X and XX 
arrays , thinning as needed . The minimum and maximum values are computed during 
the reading. Cards 92 to 96 determine the scale to be used for a plot. Cards 97 
and 98 then plot the label and axis . The time histories are plotted by cards 100 
to 105 . 

Cards 113 to 126 form the IPLT vector which lists the states, controls, and extra 
signals to be plotted . XPLOT , NUPLT , and NEXPLT define this information . The 
states , controls , and extra signals are treated together as one concatenated vector 
throughout subroutine THPLOT . States , controls , and extra signals are read into 
core and plotted up to 2 X NCH at a time by cards 128 to 168. This code is similar 
to that used for the observation plots at cards 67 to 110, with the following differences 
First, the IPLT vector is used at cards 146 and 151 to select the channels; second, 
only one line is drawn per plot instead of two; and third , the signal minima and 
maxima are initialized to 0 at cards 137 to 139 so that automatic scaling will always 
include 0. 

THPLOT is called at card MMLE3.68. 


A . 2 . 2 Utility Subroutines 

The utility subroutines perform general tasks , such as plotting and matrix 
manipulation, that are not specific to the MMLE3 program. The subroutines are 
discussed in alphabetical order. 

A . 2 . 2 . 1 ABEND . - Subroutine ABEND is an error exit subroutine. It reads 
60,000 cards from the input file in order to force an end-of-file error. This inten- 
tional error is intended to get an error traceback from the system (both IBM and CDC 
systems give traceback in response to an end-of-file error) . If an end-of-file error 
does not occur (only realistically possible if the wrong file is connected to unit 
UCARD (ref. 1, sec. 3.2)) , the subroutine stops. The return statement is never 
executed , but is required by some computers . 

A . 2 . 2 . 2 ADD (A , B , C ) . - Subroutine ADD adds the matrices A and B , placing 
the result in C . An error exit is taken if the matrix sizes are not the same . The 
physical dimension of the C matrix is assumed equal to that of the A matrix. 
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A . 2 .2 .3 ADDPAR (A , B , IBP , JBO) . - Subroutine ADDPAR adds a matrix A into 
a partition of a matrix B starting at B (IBO , JBO) . An error exit is taken if the 
partition involved exceeds the logical limits of the B matrix . 

A .2. 2. 4 AXES (XPAGE , YPAGE , AXANG, AXLEN , FIRSTV, SCALE. ANNOT , HGT) . 
Subroutine AXES draws an axis and annotation; no axis label is included . XPAGE 
and YPAGE are the coordinates of the beginning of the axis . AXANG is the angle of 
the axis from horizontal , and AXLEN is the length . FIRSTV is the beginning value 
for the annotation, and SCALE is the change per unit distance along the axis. ANNOT 
controls the placement and orientation of the annotation. Annotation can be placed 
on either the clockwise or counterclockwise side of the axis , and it can be parallel 
or perpendicular to the axis , according to the following table . 


ANNOT 

Position 

Orientation 

0 

Counterclockwise side 

Parallel 

90 

Counterclockwise side 

90° 

-90 

Clockwise side 

90° 

±180 

Clockwise side 

Parallel 

Other 

Clockwise side 

Parallel 


HGT is the height of the annotation. The variable TICDST , defined at card 13, is 
the distance between tick marks . 

Cards 23 to 42 determine the scaling exponent, NEX, and the number of digits 
placed left of the decimal point, NDECL. Cards 20 to 22 and 44 to 50 determine the 
number of di^ts to the right of the decimal point, NDECIM, and the total number of 
digits including decimal point and sign , NUMLEN . Up to five digits are allowed 
before a scaling exponent is used. Cards 52 to 70 compute starting positions and 
sign increments for the annotation . XN and YN are the starting pen positions for 
the annotation . DXNPOS and DYNPOS are increments added to the pen position if 
the annotated value is positive. Cards 72 to 83 draw the annotation values. Then 
cards 85 to 95 draw the scaling exponent if needed. Cards 97 to 110 draw the axis 
line and tick marks . Negative TICLEN puts the ticks on the side opposite the annota- 
tion . If the ticks and annotation are to be on the same side , TICLEN should be 
made positive. 

A. 2.2.5 DIGIT (I) . - Function DIGIT returns the character representation of an 
integer input argument , I . The input argument must lie in the range 0 to 20 or an 
error message will be printed and an error exit taken . The range of legal values 
can easily be modified. The ENCODE statement found on many large systems per- 
forms the task of subroutine DIGIT with much more versatility . Subroutine DIGIT 
is used, however, because ENCODE is nonstandard and, therefore, not available 
on some systems . 

A.2.2.6 DMULT (A, B, C) . - Subroutine DMULT multiplies a diagonal matrix A 
times a general matrix B and places the result in C . An error message is printed if 
the dimensions of A and B are inconsistent for multiplication or if A is not square. 
Off-diagonal elements of A are ignored if present . The physical dimension of C 
is set equal to that of B . 
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A. 2. 2. 7 EAT (A, TT, PHI, APHI, A2, A3, NEAT) . - Subroutine EAT computes 

At 

the matrix exponential e and its integral. Series expansion with time scaling is 
used (ref. 7) . Ten terms are used for the series expansion. A is the matrix which 
is to be exponentiated, TT is the time interval, and NEAT is the power of 2 used for 

A f 

time scaling . On return , PHI is the matrix exponential e , and APHI is its integral 

r* 

I 

I e ds . A2 and A3 are scratch matrices . 

•'O 


Cards 9 to 20 compute the time-scaled matrix exponential and its integral. 

Cards 22 to 32 then double the time interval NEAT times . 

A.2.2.8 EIGENG (A, Z, WR, WI, FVl, VECTS) . - Subroutine EIGENG computes 
the eigenvalues and , optionally , normalized eigenvectors of a real general matrix . 
EIGENG uses the EISPACK (ref. 6) routines BALANC , ELMHES , HQR, ELTRAN , HQR2, 
and BALBAK . These routines use the QR algorithm after balancing and elementary 
transformation to upper Hessenberg form (ref. 8) . A is the input matrix, which 
is destroyed. VECTS is an input argument, set to TRUE if eigenvectors and eigen- 
values are desired , FALSE if only eigenvalues are desired . Z is the matrix of 
eigenvalues if VECTS is TRUE , otherwise Z is not used and need not be dimensioned 
in the calling subroutine. The vectors are stored in the same order as the eigenvalues. 
For complex eigenvectors , the real part is stored in one column , followed by the 
imaginary part (corresponding to positive imaginaiy part of eigenvalue) in the next 
column . The complement eigenvector is not stored . WR and WI contain the real and 
imaginary parts of the eigenvalues. The eigenvalues are not ordered, except that 
complex conjugate pairs are adjacent, with the positive imaginary part first. FVl is 
a scratch vector the same length as WR and WI needed only if VECTS is TRUE . Error 
messages are printed if A is not square or if the QR algorithm fails . 

Cards 27 to 39 call the recommended sequence of EISPACK routines for the eigen- 
values only or eigenvalues and eigenvectors . WI is used for scratch storage and 
communication between ELMHES and ELTRAN . WR is used for scratch storage in 
BALANC; if VECTS is TRUE, WR is communicated to BALBAK via FVl. ISl and IS2 
are used for communication between the EISPACK routines. Cards 41 to 63 normalize 
the eigenvectors if eigenvectors are requested. 

A. 2. 2. 9 GET (A, MAX, II, JJ) . - Subroutine GET finds the physical and logical 
dimensions of a matrix, A. On return, MAX is the physical number of rows, II 
and JJ are the logical number of rows and columns , respectively . 

Cards 12 to 19 search for the physical dimension flag stored in the first column of 
the last row of A. In order to speed this search, a table of previously used physical 
dimensions (MAXS) is maintained. The locations of this table are checked first; 
then, locations 1 to MAXMAX are checked. Error messages are printed if the table 
of physical dimensions used becomes too long (current limit is 10) . Cards 28 and 29 
find the logical number of columns and rows . 

Subroutine GETSET must be called before the first call to GET in order to define 
TEST and MAXMAX and to initialize NUMBER to 0 . 
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A. 2. 2. 10 GETLAB (A) . - Function GETLAB returns the name of the matrix 
given as an argument , A . If the name location is 0 (common when the name was 
never defined) , the characters NONE are returned for the name. 

A ■ 2 . 2 ■ 11 GETP (A , MAX) . - Subroutine GETP finds the physical dimension 
of a matrix , A . MAX is returned as the physical dimension . If the physical 
dimension flag is not found , 0 is returned for MAX . Note that this is different from 
subroutine GET , which stops with an error message if the physical dimension flag 
is not fouhd. 

A ■ 2 ■ 2 ■ 12 GETPAR (A , B , lAO , JAO , IIB , JJB) . - Subroutine GETPAR moves 
the IIB by JJB partition of a matrix A starting at A(IA0 , JAO) to a matrix B . A and 
B may occupy the same location in storage. If the physical dimension of B was 
previously defined, it is used as defined; otherwise, the physical dimension of B 
is assumed the same as that of A . An error exit is taken if the specified partition 
exceeds the logical dimensions of A . 

A. 2.2.13 GETSET (MAX) . - Subroutine GETSET is the initialization routine for 
the matrix manipulation routines . It must be called before any of the other matrix 
routines are used. The argument, MAX, is the largest matrix physical dimension 
that will be allowed. An error exit is taken if MAX is less than 2 or greater than 101. 
GETSET defines the value , 4HTEST , used for the physical dimension pointers and 
initializes the list of physical dimensions to 0. GETSET defines the I/O unit numbers 
used by the matrix routines for the card reader (1) , card punch (2) , and line 
printer (3) , if they were not previously defined. A warning message is printed if 
GETSET defines these values . If the unit numbers are defined in GETSET and are 
inconsistent with the unit numbers that should be used, errors will generally result. 

A . 2 . 2 . 14 IDENTl (A , MAX , II) . - Subroutine IDENTl initializes a matrix to an 
identity . The first argument is the matrix , A; the second argument is the physical 
size of the matrix, MAX; the third argument is the logical number of rows, II. 

Since the matrix will be square, the number of columns is not needed. Subroutine 
ZOTl is called to set up the matrix and initialize it to 0. The diagonal elements 
are then set to 1 . 

A.2 .2.15 IDIGIT (A) . - Function IDIGIT returns the integer value corresponding 
to a single-digit, Hollerith character argument, A. A blank is treated as 0; any 
other nonnumeric character returns a value of -1. Arguments containing more than 
one digit will not be recognized; thus , they will also return a -1 value. The 
FORTRAN DECODE statement performs the function of IDIGIT with more generality 
but is not available on many systems . 

A.2.2.16 IHMSMS (ITM,T) . - Subroutine IHMSMS converts time in total 
milliseconds to hours , minutes , seconds , and milliseconds . The first argument is 
an integer with time in total milliseconds, ITM. The second argument, T, is a 
four-word integer output vector with the hours , minutes , seconds , and milliseconds . 

A . 2 . 2 .17 INV (A) . - Subroutine INV inverts a general square matrix , A , in 
place. The routine is relatively unsophisticated and will not perform well on ill- 
conditioned matrices. The algorithm used is Gauss elimination (ref. 8) with no 
pivoting. An error exit is taken if the matrix is not square. 
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A. 2. 2. 18 LINES (X, Y, NPT, ISKIP , JSKIP, HGT , L) . - Subroutine LINES 
plots solid or dashed lines through a vector of points . X and Y are the vectors of 
the X and Y values to be plotted; NPT is the number of points in these vectors . The 
locations X(NPT + 1) and Y (NPT + 1) specify the X and Y values at the origin; 
X(NPT + 2) and Y (NPT + 2) specify the X and Y scales in units per centimeter. 

ISKIP is a skipping parameter. The absolute value of ISKIP is a thinning parameter 
for the data. If ISKIP is 1, all of the data are used; if 2 , every second point is 
used, etc. If ISKIP is positive, the data are plotted starting at the beginning of 
the X and Y vectors; if negative, starting at the end. JSKIP controls dashing and 
symbols . If JSKIP is 0 , a solid line with no symbols will be drawn . If JSKIP is 
positive, a solid line will be drawn with symbols every JSKIPth point. If JSKIP is 
negative, a dashed line will be drawn with dashes JSKIP - 1 intervals long and 
spaces 1 interval long. A symbol will be put at the beginning of each dash; there- 
fore, JSKIP = -1 results in symbols only since the dashes are of length 0. HGT is 
the symbol height, and L is the CalComp symbol number (ref. 4) . 

A . 2 . 2 . 19 LYAPCB (P , A , C ) . - Subroutine LYAPCB solves a continuous-time 
block-diagonal steady-state Lyapunov equation. The form of the equation is 

AP + PA* = C 
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must be skew symmetric, and both diagonal elements of the blocks must be equal, 
must be symmetric . The first argument , P , is the symmetric solution matrix to 
the equation. The second and third arguments, A and C , are square input matrices. 
There are no checks on the structure of A and C . 


Because of the block-diagonal structure of A , the problem separates into 1 by 1 , 

1 by 2 , and 2 by 2 partitions . Each such partition is a linear equation in 1 , 2 , or 4 
elements of the P matrix. The solutions to these partitions are coded explicitly, 
taking advantage of the properties of A . 

A . 2 .2 .20 MAKE (X, Y) . - Subroutine MAKE copies a matrix. The first argument 
is the output matrix, X; the second argument is the input matrix, Y. The physical 
dimensions of the input and output matrices are assumed to be the same (see sub- 
routine MOVE otherwise) . 

A . 2 . 2 . 21 MATED (A) . - Subroutine MATED reads the body of a matrix , A , 
from the card reader file . The matrix name , number of rows , and number of 
columns must be in common block INMAT when MATED is called . The physical 
dimension of the matrix must have been defined previously . The matrix is read one 
row to a card in 8F10 format. If the number of columns is 0, the matrix is assumed to 
be diagonal , and the diagonal elements are read from one card . After reading the 
matrix , MATED calls subroutine SPIT to print it out . 

A. 2. 2. 22 MIL (T) . - Function MIL converts time in hours , minutes , seconds , and 
milliseconds to total time in milliseconds . The input argument , T , is a four-word 
integer vector . 

A. 2. 2. 23 MOVE (A,B) . - Subroutine MOVE moves a matrix A into a matrix B . 
The physical dimensions of A and B must be defined before the call and may be 
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different , MOVE differs from MAKE in the order of the arguments and the treatment 
of the physical dimensions . 

A. 2. 2. 24 MULT (A, B, C) . - Subroutine MULT multiplies a matrix A times a 
matrix B and places the result in C . The physical dimension of C is assumed equal 
to that of A . An error exit is taken if the number of columns of A is not equal to 
the number of rows of B . The CPU time spent in MULT is significant; therefore , 
it may prove worthwhile to use assembly language versions of MULT . 

t 

A . 2 . 2 .25 MULTT (A , B , C) . - Subroutine MULTT multiplies a matrix A times a 
matrix B* and places the result in C . The physical dimension of C is assumed equal 
to that of A. An error exit is taken if the number of columns of A is not equal to the 
number of columns of B . 

A ■ 2 ■ 2 .26 MVMULT (A , B , C) . - Subroutine MVMULT multiplies a matrix A 
times a vector B and places the result in matrix C . 

A. 2. 2. 27 PLOP (X) . - Subroutine PLOP punches a matrix, X, in standard matrix 
format . If the matrix name is DONT , the subroutine is bypassed . 

A ■ 2 .2 .28 PLTDAT (X , Y) . - Subroutine PLTDAT puts the date and time on a 
plot for plot identification . The arguments are the X and Y position at which the 
date and time are to be placed . Subroutine PLTDAT calls the machine-specific 
subroutines DATE and TIME to obtain the information from the system clock . 

A. 2. 2. 29 REDUCE (A, MAX, N) . - Subroutine REDUCE factors a positive definite 
symmetric matrix using Cholesky's decomposition (ref. 8) and inverts the factors. 

The arguments are the matrix. A; the physical dimension, MAX; and the logical 
dimension , N . The matrix is replaced by the factorization on exit . The matrix is 

factored into the form L ^DL* ^ , where L is lower triangular with unity diagonal 
elements and D is diagonal . On exit , the lower triangular part of the matrix contains 
L (except for the diagonal elements, which are not stored) , and the diagonal contains 
D. The strict upper triangle is not used. 

A. 2. 2. 30 RICATC (P, A, B, C, DUM, H, E, WR, WI, FVl) . - Subroutine RICATC 
solves the continuous -time steady-state matrix Riccati equation AP + PA* + B - PCP = 0. 
The first argument is the symmetric solution matrix , P . The second to fourth argu- 
ments are the square input matrices , A , B , and C . (B and C must be symmetric . ) 

The fifth argument is a dummy matrix , DUM , which is the same size as A , B , and C . 
The sixth and seventh arguments , H and E , are dummy matrices of physical dimension 
at least 2N + 1 by 2N, where N is the logical number of rows at A. The last three 
arguments, WR, WI, and FVl, are scratch 2N-vectors. 

Potter's method (ref. 9) is used for the solution. Error messages are printed 
and the program stops if the Hamiltonian has a 0 eigenvalue or if exactly half of its 
eigenvalues do not have negative real parts . 

A ■ 2 . 2 . 31 ROWCOL (IR , IC , STRING) . - Subroutine ROWCOL picks a matrix row 
and column number from a string of five characters. The first two arguments, IR 
and IC , are the integer row and column number outputs , respectively . The third 
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argument , STRING , is an input five-word vector containing one Hollerith character 
in each word. The row and column numbers may use one or two characters and 
must be separated by a comma. If less than five characters are required, any 
nonnumeric character , except a blank , can be used as a terminator . Subroutine 
ROWCOL is generally used when matrix locations are being read from cards. 

A . 2 . 2 . 32 SCALE2 (XMIN , XMAX , S , AMIN , SCALE) . - Subroutine SCALE 2 
determines reasonable plotting scales . The first two arguments , XMIN and XMAX , 
are the minimum and maximum data values , respectively . S is the axis length in 
centimeters or half inches . AMIN and SCALE are output arguments for the minimum 
value for the axis and the scale in units per centimeter or per half inch . If the data 
minimum is greater than or equal to the maximum, the scale will be set to -999, and 
the minimum value for the axis will not be defined. Scales of one, two, or five times 
a power of 10 units per centimeter or per half inch will be used by SCALE2. The 
scale and minimum value will also be chosen so that the value 0 would appear a 
multiple of 2 centimeters or half inches from the beginning of the axis (though it need 
not lie within the range of the axis) . If both positive and negative values are included , 
the axis length should be at least 4 centimeters or half inches to insure that all values 
will fit within the range of the axis . The use of centimeter or half inch units depends 
on the program's call to FACTOR (sec. 2.5) . SCALE2 is not directly affected. 

A. 2. 2. 33 SET (A, II, JJ) . - Subroutine SET defines the logical dimensions of a 
matrix. The physical dimension must have been previously defined. The three 
arguments are the matrix. A; the logical number of rows , II; and the number of 
columns , JJ . An error message is printed and an error exit taken if the logical 
number of rows is greater than or equal to the physical dimension. 

A. 2. 2. 34 SETl (A, MAX, II, JJ) . - Subroutine SETl defines the physical and 
logical dimensions of a matrix. The four arguments are the matrix. A; the physical 
number of rows , MAX; the logical number of rows , II; and the logical number of 
columns , JJ . An error message is printed and an error exit taken if the logical 
number of rows is greater than or equal to the physical dimension. 

A. 2. 2. 35 SET2 (A, MAX, II, JJ , ALAB) . - Subroutine SET2 defines the physical 
and logical dimensions of a matrix and its name. SET 2 is identical to SETl with the 
addition of a last argument , ALAB , for the matrix name . 

A. 2. 2. 36 SINV (A) . - Subroutine SINV inverts a positive definite symmetric 
matrix. A, in place. The strict upper triangle of the matrix is ignored. SINV calls 
subroutine REDUCE to compute the inverse of the Cholesky factorization (ref. 8) . 

The inverse factors are then multiplied to obtain the inverse of the full matrix . 

A . 2 . 2 . 37 SMULT (G , A , B) . - Subroutine SMULT multiplies a matrix A by a 
scalar G and places the result in B. A and B may occupy the same storage locations; 
they must always have the same physical dimensions . 

A. 2. 2. 38 SPIT (A) . - Subroutine SPIT prints a matrix, A. If the matrix name 
is DONT, the subroutine is bypassed. Matrices with more than 10 columns are 
printed in blocks of 10 or less columns at a time . 
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A. 2. 2. 39 SSIMEQ (A, B, X) . - Subroutine SSIMEQ solves the symmetric linear 
system AX = B , A is a symmetric matrix , B is a vector , and X is the solution vector . 
Subroutine REDUCE is called to obtain the inverse of the Cholesky factorization (ref. 8) 
of A . Then X is computed by multiplying B by the inverse factors. The lower triangle 
and diagonal of A contain the factorizations from subroutine REDUCE on return. The 
strict upper triangle of A is ignored . 

A . 2 . 2 ■ 40 SUB (A , B , C) . - Subroutine SUB subtracts a matrix B from a matrix A 
and places the result in C . An error exit is taken if the physical or logical dimensions 
of A and B are not the same . The physical dimension of C is assumed equal to that of 
A and B . 

A . 2 . 2 ■ 41 SUMULT (A , B , C ) . - Subroutine SUMULT replaces the matrix C by 
A*B + C . Only the lower triangular and diagonal parts of C are computed . The A 
and B matrices are assumed to have the same physical and logical dimensions; there 
is no check to verify this however . The physical dimension of C may differ from that 
of A and B . A significant amount of CPU time is spent in SUMULT , so it may prove 
worthwhile to write an assembly language version . 

A . 2 . 2 . 42 SYM (A , PRNT) . - Subroutine SYM checks a matrix, A, for symmetry. 

If the matrix is not square , an error exit is taken . If the matrix is square , but not 
symmetric, it is symmetrized using the lower triangular part. The logical variable 
PRNT controls a warning message. If PRNT is TRUE , a warning message is printed 
when A is square but not symmetric . If PRNT is FALSE , the symmetrization is done 
without comment. 

A. 2. 2. 43 SYMBL4 (X, Y, HOT, TITLE, ANGLE, NCHAR) . - Subroutine SYMBL4 
performs the same function as the standard call to the CalComp routine SYMBOL 
(sec. 2.5), except that the data are assumed to be stored only four characters per 
word. NCHAR characters of alphanumeric data, stored four characters per word in 
the array TITLE, are plotted at an angle of ANGLE degrees from the horizontal. 

NCHAR is limited to 400; otherwise, only the first 400 characters will be plotted. 

The starting position for the plot is given by X and Y . Following the usual CalComp 
convention, X and/or Y may be 999 to indicate that the plot starts at the previous 
pen position . HGT is the height of the symbols in centimeters . 

A . 2 . 2 .44 TRANSP (A , B) . - Subroutine TRANSP transposes a matrix A and 
places the result in matrix B . A and B may occupy the same storage locations . The 
physical dimensions of A and B may be different. If the physical dimension of B 
was not previously defined , it is assumed equal to that of A . 

Cards 16 to 34 handle the largest square partition of A. This part must be 
handled specially since A and B may occupy the same storage locations. Cards 37 
to 45 handle the remaining columns of A if A had more columns than rows. Cards 47 
to 54 handle the remaining rows if A had more rows than columns . 

A . 2 . 2 .45 UNSET (A) . - Subroutine UNSET deletes the physical dimension pointer 
from a matrix. A, if it is present. This has the effect of making the physical dimen- 
sion of A undefined. UNSET is used if the space in a matrix is to be used for a tem- 
porary scratch matrix with a different physical dimension. It is also used afterward, 
preparatory to redefining the correct physical dimension for the matrix . 
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A. 2, 2. 46 VMADD (A, X, B, U, S, ONES, C, V. Y) . - Subroutine VMADD forms 
the expression AX + BU + S ONES + CV and stores the result in Y . A , B , S , and C 
are matrices which may have different physical dimensions and logical number of 
columns . An error exit is taken if the logical number of rows of A , B , S , and C 
are not all the same. X, U, ONES, V, and Y are vectors. 

A . 2 .2 .47 ZMULT (A , B , C) . - Subroutine ZMULT computes the matrix expres- 
sion C = C + AB . An error exit is taken if the logical number of columns of A does 
not equal the number of rows of B . The physical dimensions of A and C are assumed 
equal. ZMULT takes advantage of O's and I's in the A matrix to save time. There- 
fore , if A is sparse , ZMULT is quite efficient . For general A matrices , subroutine 
MULT may be more efficient . A significant amount of CPU time is spent in ZMULT , 
so it may prove worthwhile to use an assembly language version . 

A . 2 . 2 . 48 ZOT (A) . - Subroutine ZOT zeros a matrix , A . The physical and 
logical dimensions of A are assumed to have been previously defined . 

A. 2. 2. 49 ZOTl (A, MAX, II, JJ) . - Subroutine ZOTl defines the physical and 
logical dimensions of a matrix. A, and then zeros the matrix. MAX is the physical 
number of rows . II and JJ are the logical number of rows and columns , respectively . 

A. 2. 2. 50 ZOT2 (A, MAX, II, JJ , ALAB) . - Subroutine ZOT 2 defines the physical 
and logical dimensions and the name of a matrix , A , and then zeros the matrix . MAX 
is the physical number of rows . II and JJ are the logical number of rows and 
columns , respectively . ALAB is the matrix name . 


A . 2 . 3 Standard Aircraft Routines 

This section describes the standard aircraft user routines . The subroutines 
are listed in alphabetical order . Each subroutine is first described in terms of its 
general function in the program. The communication to the basic program is 
described , and the place from which the subroutine is called is referenced . This 
description is intended to guide the user in coding the new set of user routines for 
a different problem . After the general description , the standard aircraft 
version of each subroutine is discussed . Subroutine WTDEF is called from MATDEF , 
and subroutines INTERP and WTTRAN are called from WTDEF . WTDEF , INTERP , and 
WTTRAN are never called directly from the basic program; therefore, they do not 
have general functional descriptions independent of the standard aircraft routines— 
they are merely parts of the standard aircraft routine's implementation of subrou- 
tine MATDEF . 

A . 2 . 3 . 1 AVERAG . - Subroutine AVERAG provides the user routines convenient 
access to the time history averages. AVERAG is called at card THDATA.lll. This 
call occurs after USERIN and before MATDEF. AVERAG is therefore convenient for 
defining variables needed for MATDEF , but not read in by USERIN . The averages 
are available in common block AVGCOM. 

The standard aircraft routine AVERAG obtains QBAR , V , ALPHA , THETA , PHI , 
and MACH if they were not read by subroutine USERIN . The values 0 and 999 , 
depending on the variable, are indications that USERIN did not read these variables. 
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The channels for ALPHA , THETA , and PHI depend on whether the case is longitudinal 
or lateral- directional . The ALPHA computation includes the corrections for upwash 
and instrument position . Cards 19 and 20 define the distances that the moment 
derivatives and instrument positions must be shifted to reference them to the flight 
center of gravity . This computation has no direct connection with the time history 
averages in the subroutine supplied; it could just as well have been done in sub- 
routine USERIN . It is envisioned , however , that the user might desire to make a 
program modification to compute the flight center of gravity as a function of some 
extra channel averages . Therefore , no computations are done with the center of 
gravity until this point so that such a modification can be readily made . 

A. 2. 3. 2 INTERP (ALPHA, NABP, ABP , lA, JA, FIA, FJA) . - Subroutine INTERP 
computes indices and factors for linear interpolation . Subroutine INTERP is a stan- 
dard aircraft routine called only from the standard aircraft routine WTDEF; therefore, 
it does not have a general function independent of the standard aircraft routines . 

ABP is a vector of length NABP containing the break point values for the inde- 
pendent variable . The values in ABP are assumed to be monotone strictly increasing; 
no check is made to verify this . ALPHA is the value of the independent variable to 
which data are to be interpolated . On return , lA and JA are the indices to be used for 
interpolation , and FIA and FJA are the interpolating factors . To interpolate a vector , 
DATA , containing dependent variable values corresponding to the independent vari- 
able values in the ABP, the expression FIA * DATA(IA) + FJA * DATA(JA) would be 
used following the call to subroutine INTERP . 

If the value of ALPHA lies outside the range of the ABP values , subroutine INTERP 
limits the output to the range given; it does not attempt to extrapolate outside of the 
range . Consequently , INTERP will work correctly when NABP is 1 . 

A. 2. 3. 3 MAKEL . - Subroutine MAKEL adds in the contribution from the L suffix 
matrices to the dimensional matrices. It is called from card DIM 1.37. On entry , 
common block DIMMAT contains the dimensional matrices except for the contributions 
from the L suffix matrices . Subroutine MAKEL should compute the L suffix matrices 
and add the contributions to the dimensional matrices . The L suffix matrices need 
not be stored separately . Common block BILIN contains the logical variable TIMVAR 
(ref. 1, sec. 3.3.8(20)) and the measured observations, controls, and extra signals. 
The L suffix matrices can be functions of the measured quantities in common block 
BILIN . The L suffix matrices should not be dependent on the unknown coefficients , 
directly or indirectly (in particular they should not depend on the computed time 
histories) . 

The standard aircraft routine MAKEL computes the matrices described in refer- 
ence 1, section 4.1.3. Cards 18 to 34 define the variables VT, ALPR, THETR, and 
PHIR . If TIMVAR is TRUE , measured values from common block BILIN are used to 
define these variables . The statement function ALPHAC corrects the measured angle 
of attack for upwash and angular rates . If TIMVAR is FALSE , values from common 
block FLCOND are used; these values can have come either from average measured 
values or from NAMELIST USER (ref. 1, sec. 4.3.3) . Cards 35 to 39 compute trig- 
onometric functions of the angles . The lateral-directional L suffix matrix terms are 
added by cards 43 to 48, and the longitudinal ones by cards 51 to 54. The variables 
DGDP and DGDT defined at cards 45 and 53 are passed through common block GRAV 
to subroutine MAKEVW . 
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A. 2. 3. 4 MAKEM . - Subroutine MAKEM computes the dimensionalization ratios , 

M suffix matrices , The values computed should be stored in common block MATRAT . 
MAKEM is called from card CALLAM.ll. On entry , the M suffix matrices in MATRAT 
will be filled with I's , Common block BILIN contains the logical variable TIMVAR 
(ref. 1, sec, 3.3.8(20)) and the measured observations, controls, and extra signals. 
The M suffix matrices can be functions of the measured quantities in common block 
BILIN , The M suffix matrices should not be dependent on the unknown coefficients , 
directly or indirectly (in particular they should not depend on the computed time 
histories) . 

The standard aircraft routine MAKEM computes the matrices described in 
reference 1 , section 4.1.3, The dynamic pressure , QT , and velocity , VT , are ob- 
tained from common block BILIN or FLCOND , depending on TIMVAR . The values 
in FLCOND can have come either from average measured values or from NAMELIST 
USER (ref. 1, sec. 4.3.3). 

A.2.3.5 MAKEVW (VB , WB , FIRST) . - Subroutine MAKEVW computes the known 
forcing functions v(t) and w(t) . It is called at cards INIT.27 and GIRL, 88, The 
formal parameter names used for v(t) and w(t) are VB and WB to avoid confusion 
with velocity and weight in the standard aircraft routines , The logical variable 
FIRST is TRUE on the first time point of each maneuver so that the subroutine can do 
any reinitialization needed . The value of FIRST should not be changed by the sub- 
routine . Subroutine MAKEVW will be called at each time point regardless of the 
variable TIMVAR (ref, 1, sec. 3.3.8(20)) . The VB and WB vectors can be functions 
of the measured quantities in common block BILIN . They can also be functions of the 
L suffix and M suffix matrices computed by user routines MAKEL and MAKEM , 
because the appropriate L suffix and M suffix matrices (time-varying or not) will 
have been defined before the call to subroutine MAKEVW . The VB and WB vectors 
should not be dependent on the unknown coefficients, directly or indirectly (in 
particular they should not depend on the computed time histories) . 

The standard aircraft routine MAKEVW computes the vectors described in 
reference 1 , section 4.1.3. The time-varying velocity from common block BILIN is 
used if TIMVAR is TRUE; otherwise, the constant average velocity from common 
block FLCOND is used . The angles of attack and sideslip used in the longitudinal 
equations are the measured values corrected for upwash and instrument position . 
Cards 37 and 64 use the variables DGDP and DGDT defined by user routine MAKEL. 
The terms on these two cards represent part of the linearizations of the gravity terms 
about the measured values . If the 0 or cp equations are not integrated (determined by 
the value of MX, the length of the state vector) the gravity terms must be evahxated 
at the measured values instead of being linearized; therefore, the terms on cards 37 
and 64 are omitted. 

A . 2 . 3 .6 MATDEF (WTFILE) . - Subroutine MATDEF defines input and matrix 
defaults. MATDEF is called at card MMLE3.62. Defaults can be defined in MATDEF 
for any of the matrices that can be read from cards . These matrices are found in 
common blocks FCOM (F matrix) , GICOM (GGI matrix) , MATIN (V suffix and APR 
prefix matrices and constraint matrix HARD) , MATRIX (N suffix matrices) , and 
SOFCOM (constraint matrix SOFT) , Data from common block AVGCOM can be used to 
define the defaults . Subroutines USERIN and AVERAG are called before MATDEF , so 
any quantities defined in these two routines can also be used (plus , of course , quanti- 
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ties defined in ONCE) . The formal parameter WTFILE is a logical variable which 
indicates whether a predicted- derivative file is available. If WTFILE is TRUE , data 
may be read from the predicted-derivative file to be used in defining the defaults . 

Subroutine MATDEF is called after the matrices have been read in from cards . 

The logical function LOADED can be used to determine whether a particular matrix 
was read from cards . The matrix is used as the argument for function LOADED , which 
returns a value of TRUE if that matrix was read from cards. Normally, the defaults 
for a matHx will be skipped if LOADED returns TRUE . MATDEF can define elements 
of a matrix, even if LOADED is TRUE for that matrix; such definitions would consti- 
tute overrides of the values read from cards . Subroutine SET can be used to define 
the logical size of a matrix. 

The constraint matrices, HARD and SOFT, require special mention. The function 
LOADED can return a value of FALSE for these matrices even if the matrices were in 
fact read from cards . In this case , the default constraints should be used in addition 
to the constraints read in . The input cards specify whether the constraints read in 
supplement the default constraints in this manner or replace the default constraints 
(ref. 1, sec. 3.3.11(6) and (7)) . In order to implement the convention described in 
reference 1, section 3.3.11(6) , the LOADED function is simply used in the normal 
manner; the default constraints are used if and only if LOADED is FALSE. The user 
should be conscious , however , that constraints read from cards can be present even 
when LOADED is FALSE . Calls to subroutine SETCON can be used to conveniently 
define the default constraints . SETCON is described in section A . 2 . 1 . 9 . 5 . 

The final values of the variables MX, MZ , MU, and MB (ref. 1, sec. 3.3.8(11) 
to (14)) in common block SIZE will not have been computed at the time MATDEF is 
called because the values may depend on the results from MATDEF. Therefore, 
these values should not be used by MATDEF unless it duplicates the logic that will 
compute the final values . The standard aircraft routine MATDEF duplicates the logic 
that determines MU and MZ at cards 27 to 31, 117, and 118. 

The standard aircraft routine MATDEF defines the defaults described in reference 1 , 
section 4.3.5. Card 33 calls WTDEF to define defaults using the predicted-derivative 
data if available. Cards 35 to 41 define the default HV for longitudinal or lateral- 
directional cases. Cards 44 to 107 define the rest of the lateral-directional defaults, 
and cards 110 to 200 define the rest of the longitudinal defaults. The GGI and F 
defaults are obtained from common block GFDEFS , previously defined by subroutine 
ONCE . The variables DCGFT , ALPHA , and V were defined by subroutine AVERAG. 

Most of the remaining variables used in MATDEF were defined by subroutine USERIN. 

The cards 116 to 119 require special mention. These cards decide whether to 
use the axis transformation to obtain the C^ derivatives in the a equation from the 

and C . derivatives in the a and a equations. If both a and a observations are 

used , the program does the axis transformation . If a is not used , the transformation 

X 

cannot be done; therefore the low a approximation is used. This approxi- 

mation is implemented by setting the variable ALPR (a used for the transformation , in 
radians) to 0 at card 119 and defining C^ = at card 132; the constraints from the 


80 



A. 2. 3, 7 


fifth row of CN , DN , and HN will automatically be ignored since the derivatives in 
these locations will not be unknown . The decision on whether to use the axis trans- 
formation or the low a approximation depends on MZ , the length of the observation 
vector . The program variable MZ will not have been finally defined when MATDEF 
is called. If the program variable MZ has the value -1, it will subsequently be 
redefined as the number of rows of the GGI matrix. Therefore, subroutine MATDEF 
must check both GGI and MZ to determine what the final value of MZ will be. The 
size of GGI is obtained by the call to subroutine GET at card 117. Note that this call 
is executed after the longitudinal default GGI matrix is defined at cards 110 and 111. 

Cards 203 to 207 override the input or default HV matrices for both longitudinal 
and lateral-directional cases. These cards delete unknown biases for observations 
that are not weighted , This forestalls the common trivial error of weighting a signal 
to 0 and forgetting to delete the signal bias from the HV matrix , The bias is , of 
course , not identifiable when the signal is not weighted . 

A. 2. 3. 7 ONCE . - Subroutine ONCE does any initialization for the user routines. 
ONCE is called only once at card MMLE3.45 at the beginning of the program, not 
for each case . Any user input which is to be read only once should be read in 
subroutine ONCE . Predicted-derivative input does not belong in subroutine ONCE; 
a separate routine, WTIN, is used for it. Subroutine ONCE is called before any of 
the other user routines , including user routine WTIN . 

The standard aircraft routine ONCE reads defaults for the GGI and F matrices. 
Cards 19 to 25 define the defaults to be used if none are read in. Card 27 reads the 
matrix header cards, and card 28 checks for the end card. Cards 30 to 45 determine 
which matrices are being read in and call subroutine MATED to read the matrix 
bodies . Common block INMAT passes the matrix name and size from the header card 
to subroutine MATED . 

A . 2 . 3 . 8 OUTPUN . - Subroutine OUTPUN punches out the estimates and related 
information as desired for derivative plotting or other programs using the estimates . 
The call to OUTPUN at card SUMOUT . 19 is controlled by the NAMEEIST variable 
PUNCH (ref. 1, sec. 3.3.8(45)). 

The standard aircraft routine OUTPUN punches the nondimensional matrices and 
Cramer-Rao bound matrices. Subroutine PEOP is called to punch the matrices. Mat- 
rices with no independent unknowns will have DONT stored for their names . Subrou- 
tine PEOP will take no action for such matrices . The standard aircraft subroutine 
OUTPUN also punches averages, standard derivations, minima, and maxima for the 
measured observations , controls , and extra signals . 

A.2.3.9 READTH (INS TAT ) . - Subroutine READTH reads the input time history 
data. READTH is called at cards THDATA.48 and 87. One time frame of data should 
be returned in common block RECRD for each call to subroutine READTH . The time 
placed in common block RECRD is a four-word integer vector in hours, minutes, 
seconds , and milliseconds . The logical variable EOFTH in common block RECRD can 
be set to TRUE to indicate that good data were not placed in common block RECRD 
because no more data were available. The variable NREC in common block INORD is 
the number of data channels to be read. NREC can be ignored if desired; it is not 
used elsewhere in the program . 
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The argument INSTAT gives information about the status of the time history 
input. Subroutine READTH should not change the value of INSTAT . INSTAT is 0 on 
the first call to READTH for a case . A test for 0 can be used to control initializa- 
tion. INSTAT is 1 when the program is searching for a start time, after the first 
call to READTH . In a multiple maneuver , there will be a search for the start time for 
each maneuver of the case after the previous maneuvers have been read. INSTAT 
is 2 when a start time has been found and data to be used are being read . 

The Variables ITM and REW in common block TAPPOS give other information on 
the input status . ITM is the total time in milliseconds of the previously read frame 
(ITM is initialized to 0 at the beginning of the program) . REW indicates when a 
rewind of the time history data file is advised. Before the first call to READTH for 
each maneuver , REW is set to TRUE if the requested start time for the maneuver is 
less than or equal to the previously read time point . On all other calls to READTH , 

REW will be FALSE . The action to be taken when REW is TRUE is up to the user and 
depends on the file structure. Usually, this will mean that the desired maneuver has 
been passed, and the time history data file should be rewound. 

The variables CARD and TAPE in common block INOPT are described in refer- 
ence 1, section 3. 3. 8 (2) . These variables are passed to subroutine READTH to 
indicate the source of time history data; they are not used elsewhere in the program . 
Subroutine READTH can obtain data from any source; the interpretation of CARD 
and TAPE is up to the user . One common usage is to have READTH create a simu- 
lated time history . 

The standard aircraft subroutine READTH reads data from cards or the time 
history input data file , UDATA , as determined by the variable CARD . The file 
UDATA can be either a tape or disk file . The file UDATA is rewound whenever REW 
is TRUE and the input is from file UDATA. The argument INSTAT is not used by 
the standard aircraft routine . 

A.2.3.10 THMOD (FIRST) . - Subroutine THMOD modifies the time history data. 
Scale factor and bias corrections can be made using the NAMELIST variables des- 
cribed in reference 1, section 3. 3. 8 (8) and (9) . More complicated data corrections 
or modifications must be done in subroutine THMOD . THMOD is called once for each 
time point at card THDATA.74, after scale factor and bias corrections to the data, 
but before any other operations . The measured time histories printed , plotted , 
or used internally in the program are the modified time histories resulting from sub- 
routine THMOD . The raw data are not retained . 

The logical argument FIRST is TRUE on the first time point of each maneuver; 
this informs THMOD of time discontinuities that may require reinitialization , The 
value of FIRST should not be changed by THMOD . The time history data are in 
common block BILIN. The modified time histories should be placed back in the same 
locations in common block BILIN . 

The standard aircraft subroutine THMOD is a null routine . 

A.2 3.11THOUT (FIRST , IT , X , Y) . - Subroutine THOUT writes the output time 
history file , UTHOUT . THOUT is called once for each time point at cards GIRL. 76 
and 120. The logical argument FIRST is TRUE on the first time point of each maneuver; 
this informs THOUT of points where initialization or reinitialization may be necessary . 
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The value of FIRST should not be changed by THOUT . IT is a four-word vector with 
time in hours , minutes , seconds , and milliseconds . X and Y are the corrected state 
estimate and predicted observation, respectively. The measured observations, con- 
trols , and extra signals are in common block BILIN . Other time history parameters 
can be obtained from common block TOGRAD if necessary . The logical vector lengths 
and physically dimensioned lengths are available in common blocks SIZE and MAXIMS . 
Subroutine THOUT can write any data desired to file UTHOUT; that file is not used 
elsewhere in the program . 

The standard aircraft routine THOUT writes the time , predicted observation , 
controls , and extra signals to file UTHOUT . The complete dimensioned lengths of 
these vectors are written . 

A, 2. 3. 12 TITPLT . - Subroutine TITPLT puts user-defined title information on 
the time history plots . TITPLT is called at card THPLOT . 44 at the beginning of the 
plot for each maneuver . The origin when TITPLT is called is 3 centimeters (or 
half inches) left of the title and at the bottom edge of the plot. The paper length can 
be obtained from common block TOPLOT. The subroutine SYMBL4 (sec. A. 2. 2. 43) 
may prove useful in subroutine TITPLT . 

The standard aircraft subroutine TITPLT is a null routine . 

A.2.3.p UINIT (X, YBIAS, UMODEL) . - Subroutine UINIT defines the initial 
condition of the state and the observation and control biases . UINIT is called at 
card INIT . 28 . The outputs of UINIT are X , YBIAS , and UMODEL . X is the initial 
state . YBIAS is a bias added to the computed observations to compare with the 
measured observations . UMODEL is the initial control of the model . The program 
will compute a control bias by subtracting the measured and model controls . The 
primary data input to UINIT is the measured data in common block BILIN . 

When UINIT is called , X is 0 , YBIAS is equal to the measured observation , and 
UMODEL is 0. If these values are unchanged, perturbation equations result. A 
common alternate choice , assuming the initial state can be suitably defined , is to set 
YBIAS to 0 and UMODEL to the measured control . 

The call to UINIT is controlled by the NAMELIST variable USERIC (ref. 1, 
sec . 3 . 3 . 8 (26) ) . If USERIC is FALSE , perturbation equations are used . 

The standard aircraft routine UINIT sets YBIAS to 0 and UMODEL to the 
measured control . The initial states are defined equal to the corresponding initial 
measured observations; the measured angle of attack and sideslip are corrected for 
upwash and vane position in order to define the initial angle of attack and sideslip 
states . The initial measured velocity is used in the position correction if TIMVAR 
(ref. 1, sec. 3.3.8(20)) is TRUE; otherwise, the average velocity is used. 

A.2.3.14 USERIN (WTFILE) . - Subroutine USERIN reads input required by the 
user routines for each case. It also does any user routine initialization for the case. 
The argument WTFILE is a logical variable that indicates whether a predicted- 
derivative file is available . The value of WTFILE should not be changed by subroutine 
USERIN. USERIN is called at card EDIT. 117 after EDIT has defined the defaults for 
NAMELIST INPUT (ref. 1, sec. 3.3.8) and the channel labels (ref. 1, sec. 3.3.9). 
Therefore subroutine USERIN can change the defaults defined by subroutine EDIT . 
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The standard aircraft routine USERIN reads the NAMELIST USER described in 
reference 1, section 4.3.3. It also changes several of the defaults set in EDIT as 
described in reference 1, section 4.3.6. 

Cards 49 and 50 redefine the USERIC and NREC defaults (ref. 1, sec. 3. 3. 8 (6) 
and (26)) of subroutine EDIT. Then cards 51 to 99 define the defaults for NAMELIST 
USER. The defaults defined at cards 71 to 90 are redefined at cards 91 to 98 from 
the predicted-derivative file if such a file is available. WTCG, defined at card 94, is 
not used 'Unless a predicted-derivative file is available. 

The NAMELIST USER is read at card 101. Card 102 forces SHIFT (ref. 1, 
sec. 4. 3. 3 (2)) to FALSE if no predicted-derivative data are available, and card 103 
forces LATR to be consistent with LONG (ref. 1, sec. 4. 3. 3(1)). Cards 104 to 106 
define the acceleration of gravity depending on METRIC and divide the weight by the 
acceleration of gravity to obtain the mass . 

Cards 107 to 154 redefine the subroutine EDIT defaults for the channel numbers 
and labels . Most of these defaults depend on whether the case is longitudinal or 
lateral-directional. These cards also process the UVAR defaults (ref. 1, 
sec. 4.3.3(22)) . Since the UVAR defaults depend on LONG, which is not known 
until after UVAR is read , the treatment of these defaults is somewhat unusual . 

UVAR was initially (cards 51 and 52) set to -999. Then cards 116, 117, 136, 153, 
and 154 apply the defaults to any UVAR elements that were not read in (those that 
are still -999) . 

Finally, cards 156 to 165 print out all of the information read by USERIN. 

A. 2.3.15 WTDEF (MUSE) . - Subroutine WTDEF obtains derivative estimates 
from the predicted- derivative file. WTDEF is called only from the standard aircraft 
routine MATDEF; therefore, it does not have a general function independent of the 
standard aircraft routines. The argument, MUSE , is the length of the control vector 
that will be used in the equations . 

Cards 18 to 24 position the predicted-derivative file and read the relevant header 
information from it. Then cards 25 to 27 compute interpolating factors and indices to 
interpolate the data to the ALPHA , MACH , and PARAM in common block FLCOND . 

Cards 28 to 32 compute DO-loop limits dependent on the interpolation indices. 

Cards 34 to 80 loop to read and interpolate the predicted-derivative tables. 

Cards 34 to 39 read a derivative header card and decide whether that derivative will 
be used. A derivative is not used if the type (LONG or LATR) is wrong or if the 
matrix where that derivative goes has been read from cards . The function LOADED , 
used in subroutine MATDEF to determine if a matrix has been read , is not convenient 
to use here since LOADED does not take the matrix name for its argument . Therefore 
cards 38 and 39 substitute for the LOADED function. Card 38 calls the function 
MATNO to obtain a matrix number from the matrix name. Then card 39 uses the matrix 
number as an index for the vector of input flags stored in common block MATLAB to 
determine if the matrix has been read . If a derivative is not used , the program goes 
to statement number 290 , which skips that derivative table and jumps back to process 
the next derivative . Cards 40 to 60 read through a predicted-derivative table and 
interpolate using the interpolating information computed earlier. Up to eight values 
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from the table are used in the interpolation (two break points each for ALPHA , MACH , 
and PARAM) , Fewer values may be used if only one break point is used for one or 
more directions , Card 42 interpolates between the two angle of attack break points 
at the first Mach number and parameter (param) break point . If two Mach number 
break points are used, card 46 interpolates between the two angle of attack break 
points at the second Mach number break point and first param break point; then 
card 47 interpolates between the two Mach number break points . If two param break 
points are used, cards 49 to 56 repeat the logic of cards 40 to 47 for the second param 
break point, and card 57 interpolates between the two param break points. Cards 58 
to 60 skip any remaining cards in the predicted- derivative table to position the file 
at the beginning of the next derivative table. Cards 61 to 77 place the interpolated 
value in the appropriate matrix location as specified on the predicted-derivative 
header card. 

After all of the values from the predicted-derivative file have been placed in the 
matrices , card 82 calls subroutine WTTRAN to do any required transformation on 
the data. 

A. 2. 3. 16 WTIN . - Subroutine WTIN reads predicted-derivative data from cards 
and writes a predicted- derivative file. WTIN is called at card MMLE3.53 after user 
routine ONCE but before any of the other user routines. 

The standard aircraft subroutine WTIN reads the data in the form described in 
reference 1, section 4.3.2 and writes the predicted-derivative file in the form 
described in reference 1, section 4.2.2. The primary processing involved is to 
reorder and expand the simplified data tables read in from cards . The expanded and 
reordered forms are simplest to interpolate for later use , but it is preferable to 
allow more flexibility in the input formats . 

Cards 27 to 30 read and write the title card for the data set. Cards 19 to 21, 

31, and 32 define defaults for the variables in NAMELIST WIND, and card 33 reads 
the NAMELIST. These data are written on the predicted-derivative file by cards 38 
to 42 . The break points are read from cards and written on the predicted-derivative 
file by cards 44 to 49 . Optional printed output is done by cards 50 to 61. 

Cards 63 to 126 constitute the primary operations of the standard aircraft routine 
WTIN . These cards read a derivative data table from cards , reorder and expand the 
table in triply dimensioned array BDAT , and write the data to the predicted- 
derivative file and the line printer file . 

This part of the subroutine repeats for each derivative table until an end card is 
found. Cards 63 and 64 read a derivative header card and test for the end card. 
Subroutine ROWCOL is called at card 65 to translate the character string "SUB" into 
integer row and column numbers . This information was read in as a character string 
to allow more freedom in the input format than the FORTRAN I format specification . 
Subroutine ROWCOL is described in section A.2 .2 . 

The derivative functional dependence , FS , is described in reference 1 , 
section 4 . 3 . 2 . 6 . Cards 66 to 69 set FS to the default , AMP , if blank was read in . 

The variables II, Jl, Kl, 12, J2, and K2 are used to expand the derivative table. 

II to 12 are the indices of the angle of attack break points to which a single input 
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value will be copied. If the input data are not a function of angle of attack (i.e. , 
none of the FS values are A) , II will be 1 and 12 will be NABP so that each data 
point read will be copied to all of the angle of attack break points . If the input data 
are a function of angle of attack, II and 12 will both start at 1, indicating that the 
first data point read will be copied only to the first angle of attack break point; for 
subsequent data points , II and 12 will be incremented whenever a new angle of attack 
break point is being processed. Similarly , J1 and J2 relate to Mach number break 
points , and K1 and K2 relate to param break points . 

Incrementing is done in the order specified by FS . First , the dimension 
specified by FS(1) is incremented, corresponding to different fields on one input 
card . After all of the break points of the FS (1) dimension are done , the FS (2) di- 
mension is incremented by 1, and the FS(1) break points are redone for the second 
FS (2) break point. FS(3) is incremented as the outermost loop. The variable 
INDEX keeps track of which dimension is being incremented . When all of the data 
have been read in , the expanded array BDAT is written on the predicted- derivative 
file and optionally printed. The program then loops back to process the next deriva- 
tive table . 

A.2.3.17 WTTRAN (AXIS , MUSE) . - Subroutine WTTRAN does various trans- 
formations on the data obtained from the predicted-derivative file . WTTRAN is 
called only from the standard aircraft routine WTDEF; therefore, it does not have a 
general function independent of the standard aircraft routines . 

The standard aircraft routine WTTRAN transforms longitudinal stability axis 

derivatives to body axes, computes C., and C . from total C.. and C . , and trans- 

Nq Aq n a 

forms moment derivatives from reference to flight center of gravity . The argument 
AXIS should be either "STAB" or "BODY" to indicate the axis system (stability or 
body, respectively) of the longitudinal data. The argument MUSE is the number of 
controls to use in the transformation . If the longitudinal data are already in body 
axes , the transformation from stability to body axes is skipped . The lateral- 
directional data are always assumed to be in body axes , regardless of AXIS . The 
logical variable SHIFT in common block INERTS controls the center of gravity 
transformation . If SHIFT is FALSE , the moment derivatives are assumed to be already 
referenced to the actual flight center of gravity , so no transformation of them is done; 
any difference between the values given for the flight and reference center of gravity 
is ignored. 

Just as in subroutine MATDEF, the LOADED function is used to determine if a 
matrix has been read in from cards . Subroutine WTTRAN does not change any 
matrix that was read from cards; its transformations are intended only for defining 
default matrices using the predicted-derivative file . 

Cards 20 to 36 do the transformation from stability to body axes for the longitudinal 
data . The transformation for the angle of attack derivatives on cards 34 to 36 assumes 
that total and are in HN (4,1) and HN (5,1). Therefore , this code must be after 

the transformation of and to and (cards 29 to 32) and before the compu- 
tation of and C . to replace C., and C . (cards 39 to 43) , 
iVfl Aq n a 
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Cards 39 to 43 compute the extended C.. and C . from total C., and C . by 

JVq Aq n a 

subtracting the linearized contributions of the angle of attack and control derivatives . 
The pitch rate contribution is not subtracted because the total and predictions 

are assumed to be given for zero pitch rate , even though the maneuver may actually 
have a significant average pitch rate. Note that the Cj^ and C . computed are 

functions of the angle of attack. They are consistent with the usual definition of 
Cj^ and C . only if a is 0. We often refer to the program's quantities as and 
0 0 0 
C. extended, to distinguish them from the usual definition. 

^0 


Center of gravity transformation is done by cards 48 to 55 for lateral-directional 
data, and cards 58 to 65 for longitudinal data. The longitudinal transformations 
assume that the derivatives are in body axes; therefore , this code must follow the 
stability-to-body axis transformation. 


A . 2 . 4 EISPACK Routines 

The EISPACK routines used by MMLE3 are BALANC, BALBAK, ELMHES , ELTRAN, 
HQR , and HQR2 . These subroutines are exactly as obtained from Argonne Labora- 
tories and described in reference 6, with one exception. Cards 103 and 104 of HQR 
and cards 127 and 128 of HQR2 have been modified to increase the maximum number 
of QR iterations from 30 to 50. We have found cases where the routines require 
more than 30 QR iterations to converge to the specified accuracy on CDC computers . 
The individual subroutines are not described here , as adequate documentation is 
provided by reference 6 . The EISPACK routines are called only from subroutine 
EIGENG. 

The variables RADIX and MACHEP in subroutines BALANC, HQR, and HQR2 
are machine dependent. These variables are discussed in section 2.4. 
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APPENDIX B 
PROGRAM COMSUB 


Program COMSUB performs common deck substitution , Users with access to 
CDC UPDATE (ref. 2) or other similar utility packages will not need to use COMSUB. 

It reads the UPDATE source cards described in section 1.1, substitutes the common 
decks in appropriate places , and punches the resulting FORTRAN decks . 

No input cards are necessary except for the UPDATE source cards . 

The program card and cards 43 to 46 of COMSUB can be altered as desired to 
define the file numbers for the card reader (UCARD) , card punch (UPUNCH) , line 
printer (UPRINT) , and UPDATE source cards (UDATA) . COMSUB does not use the 
file UCARD . The file UDATA can be assigned to the card reader , a disk , or tape 
file , depending on how the source data are available . 

The maximum number of common decks allowed is 100 , and no more than 200 total 
cards are allowed for all of the common decks . These dimension limits are easily 
charged on the dimension statement at card 38 and in the definitions at cards 47 and 48 . 

Program COMSUB requires an end-of-file check . The comment cards describe 
how to configure the program for IBM or CDC end-of-file conventions . Only three 
cards are affected . Alternately , the program can omit the end-of-file checks and 
watch for a card with ”*END" in the first four columns to flag the end of the source 
deck . This alternative is not machine dependent , but does require that the *END 
card be at the end of the source data . 

Cards 54 to 60 initialize the common deck counts to 0 and read the first card . 

This first card must be a *COMDECK card , or the program will do an error stop . 

Cards 62 to 69 start a new common deck . The number of common decks is incre- 
mented by 1 and checked against the maximum . The common deck name is stored in 
CNAMES . A pointer in ICARDl is defined to point to the first column in the matrix 
CDECKS included in this common deck. 

Cards 71 to 79 read a common deck, storing each card in a column of the matrix 
CDECKS . A dimension limit check is made on the number of columns allowed . Any 
card with a star in column 1 will be assumed to be an UPDATE directive , and thus 
will define the end of the common deck . 

Cards 81 to 87 finish the storing of a common deck and decide whether more 
common decks follow . A pointer in ICARD2 is defined to point to the last column 
in CDECKS included in the common deck. The UPDATE directive that follows the 
common deck is then examined. If it is a *COMDECK directive, another common deck 
follows , so the program loops back to read the new common deck . If it is a *DECK 
directive , the common decks have all been read , and the program continues to 
the next section . Any directive other than *DECK or *COMDECK will be flagged as 
an error and the program will stop . 
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Cards 89 to 117 read and copy the remaining source cards, substituting common 
decks as called. Card 91 is used in the IBM version; card 92 is used in the CDC and 
machine general versions to read a source card. Card 93 is the CDC end-of-file 
check. When IBM end-of-file is detected, the program jumps to statement 800, and 
a card still remains to be processed; therefore, the variable LAST is defined at 
card 119 , and the program jumps back to finish processing the last card. 

If the card that is read does not have a star in column 1, cards 95 to 97 copy 
it to the punch file and loop back to read the next card. Cards 94 and 98 to 101 
decide on the processing of source cards beginning with a star . If an *END card 
is found , processing is terminated . If a *DECK card is found , the deck name is 
printed out , and the card is otherwise ignored . If a *CALL card is found , common 
deck substitution is done. Any other card beginning with a star is copied to the 
punch file without special treatment . 

Cards 103 to 117 perform common deck substitution when a *CALL card is 
detected . Cards 104 to 107 search the list of common deck names for the requested 
common deck . If a matching name is not found , the program stops with an error 
message. Cards 111 to 117 punch a copy of the common deck and jump back for 
processing of the next source card. Note that zero length common decks are allowed, 
in which case no cards are punched for them . 
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PROGRAM COMSUB(lHPUT-258,PUNCH-512fOUTPUT«258fOATA-512# 

CGMSU6 

2 


TAPEl- INPUT# TAPE2-PUNCH,TAPE3-OUTPUT,TAPEA-DATA) 

COhSUB 

3 

c 


COMSUB 

A 

c 

RICHARD E. PAINE 9 JAN 79 

COMSUB 

3 

c 

PROGRAM TO CREATE FCRIRAN SOURCE DECKS FROM UPDATE SOURCE DECKS, 

COMSUB 

6 

c 

READS IN COMMON DECKS FOLLOWED BY REGULAR DECKS, 

COMSUB 

7 

c 

SUBSTITUTES COMMON DECKS INTO REGULAR DECKS AS CALLED FCR, 

COMSUB 

8 

c 

PUNCHES OUT RESULTING SUBSTITUTED DECKS. 

COMSUB 

9 

c 


COMSUB 

10 

c 

END-OF-PILE IS CHECKED BY ONE OF THREE METHODS, 

COMSUB 

11 

c 

1) CDC TYPE EOF FUNCTION (INDICATES KQ DATA WAS RETURNED ON 

COMSUB 

12 

c 

PREVIOUS READ) 

COMSUB 

13 

c 

FOR THIS TYPE CFtCK# THE CARO FLAGGED IBM IN COL, 73 SHOULD 

COMSUB 

14 

c 

HAVE A C IN CCL. 1 SINCE COC COMPILERS CO NOT RECCGNIZE 

COMSUB 

13 

c 

THE SYNTAX CF THIS CARD, THE CAROS FLAGGED WITH CDC IN 

COMSUB 

18 

c 

COL, 73 SHOLLD LEAVE CLL* 1 BLANK. 

COMSUB 

17 

c 

2) IBM type end- parameter (BRANCHES IF DATA JUST RETURNED WAS 

COMSUB 

18 

c 

LAST RECORD IN FILE) 

COMSUB 

10 

c 

FOR THIS TYPE CHECK# THE CARDS FLAGGED CDC IN COL, 73 SHOULD 

COMSUB 

20 

c 

HAVE A C IN CCL. I AND THE CARO FLAGGED WITH IBM IN COL. 73 

COMSUB 

21 

c 

SHOULD Be BLANK IN COL* 1, 

COMSUB 

22 

c 

3) FOR ANY MACHINE - PROGRAM vHECKS FOR A DATA CARO CONTAINING 

COMSUB 

23 

c 

♦END STARTING IN COL. 1, 

COMSUB 

24 

c 

THIS CHECK IS ALWAYS VALID# BUT REQLIRES THE APPROPRIATE 

COMSUB 

25 

c 

CARL TO BE AT THE END OF THE DATA DECK. THc PROGRAM CHECKS 

COMSUB 

26 

c 

FOR THIS CONDITION IN EITHER THE COC OR IBM CONFIGURATIONS, 

COMSUB 

27 

c 

FUP A MACHINE GENERAL PROGRAM THAT USES ONLY THIS CHECK# 

COMSUB 

26 

c 

PUT A C IN COL, 1 OF CARDS FLAGGED WITH EITHER COC CR IBM 

COMSUB 

20 

c 

IN COL* 73 XCEPT FOR THE ONE CARO ALSO FLAGGED WITH ANY, 

COMSUB 

30 

c 


COMSUB 

31 

c 

LIMITATIONS# DIMENSION LIMITS 

COMSUB 

32 

c 

MAX NUMBER OF C0PM(jN DECKS IS IOC (NCOMMY) 

COMSUB 

33 

c 

Max total number Cf Cards in common decks is ZuO (NCROMX) 

COMSUB 

34 

c 


COMSUB 

35 


INTEGER UCARO#UPUNCH#UPRINT#UOATA 

COMSUB 

36 


LOGICAL LIST 

COMSUB 

37 


DIMcNSlON C&ECKS(2I#2C0 )#CNaMES(2#IOO)#ICARC1(1OO)#ICARD2(1OO)# 

COMSUB 

38 


CaMCR0(21)#CAR0(22) 

COMSUB 

30 


DATA STAR/1H»/#C0MD#ECK/AHC0M0#3HECK/#LECK/AHDFCK/#CALL/AHCALL/# 

COMSUB 

40 


EN0/3HEND/ 

COMSUB 

41 

c 


COMSUB 

42 


UCARD • 1 

COMSUB 

43 


UPUNCH • 2 

COMSUB 

44 


UPRINT - 3 

CCMSUB 

45 


UOATA - A 

COMSUB 

46 


NCROMX • ?ac 

COHSUB 

47 


NCUMMX « 100 

COMSUB 

48 


NCPUN • 0 

COMSUB 

40 


Last • .FALSE, 

COMSUB 

50 


WRiTF(UPRlNT#3000) 

COMSUB 

51 

c 

read COMMCN DECKS, 

COMSUB 

52 

c 

♦*♦♦♦♦*♦♦♦ INITIALIZE AND READ FIRST CARO. 

COHSUB 

53 


NCOhS • 0 

CCMSUB 

34 


NCARO - 0 

COMSUB 

33 


REa0(U0ATA»100I) comcrd 

COMSUB 

56 


IF(C0MCRD(1)*E0.STAR .AND* 

COMSUB 

57 


CaMCR0(2),EQ,C0HD .AND, COMCRO ( 3 ) .EQ.ECK ) GO TO 100 

COMSUB 

58 
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URlTEiUPI^lNTf 30Cb) COPCRO 
STOP 

C STAR! NEh CCPHCjN DECK, 

ICO NCOHS • NCO^S^l 

ICAkOKNCCHS) • NC*ROtl 

WRlTECUPRIHTf 30C1 ) hCCPS*CCKCRD ( A ) t CUHCRO ( 5 ) # IC ARDl (NCOHS ) 
IF(NCC«S.LE.NCOPHX) 6C TC 110 
WKlTE{UPPlNT#3wC2) NCCPhX 
STCP 

110 CNAr^ES(I#NCChS) - CCMCRD(A) 
tNAME5(3#>»CC«S) • CCHCR0(5> 

C RE40 TO ENC CF COflHON DECK, 

2C0 READ(U0ATA#100n CQPCRD 

IF(CaPCKOCl),EO«$TAR) CL TO 3wO 
NCaRO • NCAkO^l 
IF(NLARD«LE,NCRLHX) 6L 10 21u 
i»RITE(UPRINT#30C3) NCROPX 
STOP 

210 00 220 1-1*21 

22n COECKS(I#NCARO) • CCHCRO<I) 

GO TO 2uO 

C EHO CF CCHPLN DECK, 

300 ICAR02(NC0HS) - NCARO 

IF(CaPCR0(2).cQ,LCHC •AND, COHCRO< 3 ) • EQ,E CK ) GO TC 100 
IF(COPCRO(2).cQ,DECK) €C TO 350 
URIIE(UPPlNT*30t^) CGPCRO 
STOP 

350 CARD(4) - COHCRCO) 

IAR0(5) - CCiHCRC(A) 

c pead and copy oecks* 

400 URlTE(UPRINT*3vC5) CARD (4) *CARO( 5) 

420 IF(LAST) go to 90C 

C R£AC(U0ATA*iaC2#END-6OO) CARD 

READ(U0ATa* 1OO2) CARO 
IFCEaFfUDATA),NE.O.) GO TO 900 
430 IF(CARD(n*EQ,STAR) 6C TO 450 
440 NCPUN • NCPl'N+1 

URXTE(0PUNCH#1UL2) CARP 
60 TO 420 

450 IF(CARD(2),E0,ENDI GO TC 9w0 
lF(CARC(2)«tQ«0ECKI 60 TO 400 
IF(CaR0(2),E0,CALU 6L TO 5U0 
GO TO 440 

C ♦*♦♦♦♦♦♦♦♦ SEARCH FCR CALLED COPMON DECK, 

500 WR1TEUJPR1NT>30C6) CARC (4 ) ^CaRD (5 ) 

00 550 1C0M-1#NC0HS 

IF(CARD(4),EQ,CNAPE5(1*1CQH) ,ANO« CARO (5 ) ,EQ,CNAHES( 2* ICCH ) I 
GO TC 600 
550 CONTINUE 

WRITE(UPRINT>3u07) 

STOP 

C ♦♦♦♦♦♦♦♦♦♦ PUNCH CALLED COHHUN DECK, 

600 NCCROS • lCAR02acQh)*lCAROl(iCON)4l 
1F(NCCR0S*EQ,U) GL TO 420 
NCPUN - NCPUN4NCCR0S 
ICROI - ICAROKICCH) 

ICR02 • ICAR02(1CGN) 



CQPSUB 

59 


COPSUB 

60 


CCPSU6 

61 


CQPSUB 

62 


CQPSUB 

63 


CCPSUB 

64 


CQPSUB 

65 


CCPSUB 

66 


CQPSUB 

67 


CQPSUB 

68 


CQPSUB 

69 


CCPSUB 

70 


CCPSUB 

71 


CQPSUB 

72 


CCPSUB 

73 


CQPSUB 

74 


CQPSUB 

75 


CCPSUB 

76 


CQPSUB 

77 


CQPSUB 

78 


CQPSUB 

79 


CQPSUB 

80 


CCPSUB 

81 


CQPSUB 

62 


CQPSUB 

83 


CCPSUB 

84 


CCPSUB 

85 


CQPSUB 

66 


CCPSUB 

87 


CQPSUB 

68 


CQPSUB 

69 


CQPSUB 

90 

18N 

CQPSUB 

91 

CDC 

ANY CQPSUB 

92 

COC 

CQPSUB 

93 


CQPSUB 

94 


CQPSUB 

95 


CQPSUB 

96 


CCPSUB 

97 


CQPSUB 

98 


CQPSUB 

99 


CQPSUB 

100 


CQPSUB 

101 


CQPSUB 

102 


CQPSUB 

103 


CQPSUB 

104 


CCPSUB 

105 


CCPSUB 

106 


CQPSUB 

107 


CQPSUB 

106 


CQPSUB 

109 


CCPSUB 

lie 


CQPSUB 

111 


CQPSUB 

112 


CQPSUB 

113 


CQPSUB 

114 


CQPSUB 

115 
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WPITcClPUNCH^^OCl ) ((CCECKS( J»l)#J-l/21)#l-lCRDl«lCf«D2) 

GO TO ^20 

C END-OF-FILE PROCESSING. 

000 LAST - .TRUE. 

GO TO A30 

9C0 WRITECUPRINTySOlO) hCPUN 
C 

1001 F0RHAT(Alfl9A4» A3) 

1002 FORNATC Al#AA«Al#16AA#A2) 

3000 FCRMATCICOHSUB PROGRAN FOR INSERTING CCNMON DECKS"/) 

3001 FORMAK" CQNNON OECK">I AUX, 2AA# " STARTS AT CARD"# I5» ". " ) 
j002 FORNATC"OR*R ERROR. LIMIT 0F"#I5#" COMMON DECKS REACHED.") 

3003 FURMAT("O^RRP ERROR. LIMIT 0F"#16»" TOTAL CAROS IN COMMCN DECKS "» 
"REACHED.") 

30CA F0RMAT("0PAAA error. DIRECTIVE NOT RECOGNIZED - "# AI# 1<>AA#A3) 

3005 FORMAT("OCECK "#2AA»" BEING PROCESSED.") 

3006 FORMAT(" CALLING FOR tCMMON DECK "#2AA) 

1007 FURMAT("OPP* ERROR. ABOVE CQMMUN DECK NOT RECOGNIZED.") 

10C8 FLRHATt"OP«^ ERROR. FIRST CARD IS NOT A VALID COMDECK CARD."/ 
iX#Al#19A4#A2) 

3010 FQRMAT«"0T0TAL CF"#I6#" CARDS PUNCHED.") 

STOP 

END 


COMSUM 116 
COMSUB 117 
CQMSU6 118 
COMSUB 119 
COMSUB 120 
CCMSUB 121 
COMSUB 122 
COMSUB 123 
COMSUB 12A 
COMSUB 125 
COMSUB 126 
COMSUB 127 
CCMSUB 128 
COMSUB 129 
COMSUB 130 
COMSUM 131 
CCMSUB 132 
CCMSUB 133 
COMSUB 134 
COMSUB 135 
COMSUB 136 
COMSUB 137 
COMSUB 138 


c 


APPENDIX C 
PROGRAM COMPUN 


Program COMPUN punches common decks for the MMLE3 program, substituting 
in the desired physical dimensions. Used in conjunction with program COMSUB 
(appendix B) or CDC UPDATE (ref. 2) , COMPUN allows the physical dimensions of 
the MMLE3 program to be easily changed as discussed in section 3.2. 

The file numbers used by COMPUN can be altered by changing the program card 
and cards 61 to 64. The input to COMPUN is on two separate files. The card reader 
file (UCARD) contains the NAMELIST IN, which defines the physical dimensions to 
be used . The variables in NAMELIST IN are described in the comment cards . 

The file UDATA contains a template for the common decks , UDATA can be 
assigned to the same file as UCARD , in which case the template would follow the 
NAMELIST in the input cards . The template consists of the common decks (complete 
with *COMDECK cards) , except that symbols are used for the dimensions. COMPUN 
will copy the template to the punch file , with the appropriate values substituted for 
the symbols. 

Each symbol used in the template consists of a star followed by a two-character 
name. The two-character names used are described in the comment cards. COMPUN 
will substitute a two-digit value for the two-character name and a blank for the star . 

If more than two digits may be required for the value , a star should be placed after 
the two-character name in the template as well as before it. A four-digit value will 
then be substituted in the locations of the two-character name and the two stars . The 
template must have no stars after column 8 except for those used as described above . 
The end of the template is indicated by a card with "*END" starting in column 1. 

Cards 66 to 79 define the default dimensions and then read and write the NAME- 
LIST . The default values for the dimensions are the values used with the MMLE3 
program as supplied. Cards 81 to 87 check several limitations on the allowable dimen- 
sions. These limitations are described in the comment cards. Cards 89 to 94 define 
dimensions computed from the basic dimensions read in NAMELIST IN . The loop at 
cards 97 to 107 converts each of the dimension values to four Hollerith digits. 

Cards 110 and 111 read a card of the template and check to see if it is the *END 
card. Cards 120 to 134 insert two- or four-digit values in the places specified by 
the template . 

The program listing for COMPUN is shown below , followed by a listing of the 
template . 
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oooooor>or>oor>ooo«nrtr>r>or>or>oor»or)r>or>oooorkoo^o 


PROGRArt CCMPUNCCCJPPUN#RE$nE» INPUT# OUTPUT# 


COMPUN 

2 

TAPE4-CaNPUN>TAPE5«RESIZE# 


CQMPUN 

3 

TAPE 1-INPUT# TAPE 3-LlTPUT) 


COMPUN 

A 



CQMPUN 

5 

PUNCHES HNlE3 CUWMCN DECKS OR UPDATES 


COMPUN 

6 

PkOGRAH LAST MODIFIED 15 AUG 79 RICHARD HAIHE 


COMPUN 

7 



COMPUN 

8 

INPUT VARIABLES IN NAPELiST /IN/ 


COMPUN 

9 

CCORReSPONDING Z CHARACIER NAMES IN PARENS)* 


COMPUN 

XO 

MAXX - MAX NC CF STATES 

(X ) 

COMPUN 

XI 

HAXZ • MAX NL CF OBSERVATIONS 

(Z ) 

CCMPUN 

X2 

MAXU - MAX NL CF CCNTRUL INPUTS 

(U ) 

COMPUN 

X3 

HAXB - MAX NC CF BIAS INPUTS 

(B ) 

COMPUN 

XA 

LEX - MAX NO OF EXTRA SIGNALS 

(EX) 

CCMPUN 

X5 

NT - MAX NC OF INCEPENOENT UNKNOWNS ^ Z 

(HI) 

CCMPUN 

X6 

MAXTV - MAX NO OF UNKNOWNS AND CONSTRAINTS 

(TV) 

COMPUN 

X7 

MAXKV - MAX NO QP UNKNOWNS AFFECTING K 

(KV) 

COMPUN 

X8 

MAXHRO - MAX NC OF HARO CONSTRAINTS 4 1 

(HC) 

COMPUN 

X9 

MAXSFT - MAX NC CF SOFT CONSTRAINTS ♦ 1 

(SF) 

CCMPUN 

20 

NTPLT - MAX TIME POINTS FOR PLOTTING 4 2 

(TP) 

COMPUN 

21 

NSP - N'J OF SIGNALS IN CORE FOR PLrTTING 

(S ) 

COMPUN 

22 



CCMPUN 

23 

COMPUTED VARIABLES AND 2 CHARACTER NAMESi 


COMPUN 

2A 

(XI) • MAXX41 


CCMPUN 

29 

(Zl) • MAXZ41 


COMPUN 

26 

(UX) • MAXU4LEX 


COMPUN 

27 

(PL) ■ MAXX4MAXU4LEX 


COMPUN 

28 

(SI) • NSP41 


COMPUN 

29 

(S2) ■ 2*NSP 


COMPUN 

30 



COMPUN 

3X 

LIMITATIONS 


COMPUN 

32 

BECAUSE OF PROGRAM ASSUMPTIONS IN USING SCRATCH MATRICES. 

COMPUN 

33 

MAXZ «GE. MAXX 


CCMPUN 

3A 

BECAUSE OF MINIMUM SPACE REQUIRED FOR MATRIX STORAGE 

CONVENTIONS. 

COMPUN 

39 

MAXX .GE« A 


COMPUN 

36 

maxz «ge. a 


COMPUN 

37 

MAXU «GE. a 


COMPUN 

38 

MAXB .r-E. A 


COMPUN 

39 

NI .C-E. A 


CCMPUN 

AO 

FuR STATE NOISE CASE CNLY# IN SUBROUTINE KALhAHj 


COMPUN 

AX 

COHPUH prints a WARNING FOR VIOLATING THIS# BUT DOES 

NOT STCP. 

COMPUN 

A2 

HI .GE. AAMAXX4A 


CCMPUN 

A3 



COMPUN 

AA 

INTEGER UCARO#UPP1NT#UDATa#UPUNCH 


CCMPUN 

A9 

INTEGER A(66)#B(2#2C)#V(A#20)#VAL(20)#NUM(10) 


CCMPUN 

A6 

INTEGER BLANK#END#$TAR 


COMPUN 

A7 

EQUIVALENCE ( MA XX# V AL ( I ) ) # ( MAX Z # VAC ( 2 ) ) » ( M AXU# VAC ( 3 ) ) 

f 

CCMPUN 

A6 

(MAXB# VAK A) )# (LEX#VAL(5) )# (hI»VAL (6 ) ) # (MAXTV# VAC ( 7 ))# 

COMPUN 

A9 

(MAXKV# VAL (6))#(MAXHRD# VAU9))#(MAXSFT#VAL(10))» 


COMPUN 

90 

(NTPLT#VaL (^Xn# <NSP#VAL(12)) 


COMPUN 

51 

DATA STAR/IHA/# BLANK/IH /# END/AH4END/# 


CCMPUN 

52 

NUM/lHO# XHl#iH2#lH3#XHA#lH5#lH6#XH7#lHe#lH9/# 


CCMPUN 

53 

B/1HX#1H #1HZ#1H #1HU#1H #1HB#1H # IH E» 1HX#XHN# IHI# IHT# IHV# 

CCMPUN 

5A 

IHK#1HV#1HH#1H0#1HS#IHF#1HT#1HP#1HS#1H #lHX# XHX# XH Z# IHl# 

CCMPUN 

59 

XHU#lHX#iHP#lHL#lHS#XHl#XHS#XH2#AAlHt / 


CCMPUN 

56 

NAMELIST /IN/ M AXX# MA X Z# MAXU# MAXB# L E X# NI# MAXTV# HA XK V# 

MAXHRD# 

CCMPUN 

57 

MAXSFT#NTPLT#NSP 


CCMPUN 

58 


97 


c 


c 

Q FILE NUHeERS* 

UCAKO • 1 
UPRINT • 3 
UOATA ■ A 
UPUNCH • 5 

C defaolts and input* 

MAXX ■ 7 
HAXZ ■ 8 
PAXU • 4 
MAXB • 4 
LEX • 4iO 
NI . 35 
HAXTV • 50 
HAXKV ■ 15 
HAXHRO • 36 
HAX6FT • U 
NTPLT • 120? 

NSP • 3 

REAO(UCAkn»lN) 

WRITE(UPkIHT#IN) 

C CHECK LIMITATIONS* 

IF(MAXX.GT*PAXZ ) GO TC 2tO 
IF(MAXX.LT.4) GC to 200 
IF(MAXT*LT*4) Gl TO 2C0 
IF(MAXU.LT*4) GC TO 200 
IF(HAXB*LT*4) GC TO 200 
IP(NI *LT*4) Gc to 200 
IF(NI*LT.4AMAXX+4 ) WRITE (UPRINT* 3004) 

C DEFINE COMPUTED QUANTITIES. 

VAL(13) - MAXX41 
VaI(14) • HAXZ^l 
VALU5) • MAXU+LEX 
VAL(16) • MAXX^MAXU4LEX 
VAL(17) • NSPfl 
VAK16) • 2«NSP 
N6 - 16 

c CONVERT TO CHARACTER. 

DC 10 I-1*NB 
U-VAL(I)/1U0U 
J»VAL (D-Tl^iOOO 
I2»J/100 
J«J-I2»100 
T3-J/10 
I4«J-I3^10 
V(1*T )«NUM(Ii + l ) 

V(2*I )«NUMn?*l) 

V(3*I )«NUM(I3 + 1 ) 

10 V(4#i)«NUM(l4^1) 

C 

C ^t^^^tt****t*^*^***4*t4^^^4*** RfcAO TEMPLATE. 

100 READ(UOATA*1001) A 

1F(A(1)*EQ*EN0) GO TC 500 
I • 2 
110 I-I41 

1F(1*LT«64) GO TO 120 
WRITE(UPUNCH*1001) A 


CCMPUN 59 
COMPUN 60 
CCPPUN 61 
COMPUN 62 
COMPUN 63 
CCMPUN 64 
COMPUN 65 
COMPUN 66 
CCMPUN 67 
COMPUN 66 
COMPUN 69 
COMPUN 70 
COMPUN 71 
CCMPUN 72 
COMPUN 73 
CCMPUN 74 
COMPUN 75 
COMPUN 76 
COMPUN 77 
COMPUN 76 
COMPUN 79 
COMPUN 60 
COMPUN 81 
COMPUN 82 
COMPUN 83 
CCMPUN 84 
CCMPUN 85 
COMPUN 66 
CCMPUN 87 
COMPUN 88 
CCMPUN 89 
COMPUN 90 
COMPUN 91 
COMPUN 92 
COMPUN 93 
COMPUN 94 
COMPUN 95 
COMPUN 96 
COMPUN 97 
COMPUN 98 
COMPUN 99 
CCMPUN 100 
COMPUN 101 
CCMPUN 102 
COMPUN 103 
CCMPUN 104 
CCMPUN 105 
CCMPUN 106 
COMPUN 107 
COMPUN 106 
COMPUN 109 
CCMPUN lie 
COMPUN 111 
COMPUN 112 
COMPUN 113 
COMPUN 114 
CCMPUN 115 
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c 


WRITE<UPRlNT>a002) A 
GO TO lOw 

12C lF(A(n.N£.STAR) GO TC AXO 

C I^SERT QIHENSTONS WHERE INDICATED. 

DO 130 

IF(A(I*1).NE.6(1»J) .CR. A(1^2).NE«R(?#J)) GO TO 130 

IF (A(I>3).£C.STAR ) GO TC 125 

Aa)>BlANK 

A(I41)*V(3#J) 

A(I^2)«V<A»J) 

I-I + 3 
GO TO 110 
125 A(n*va>j) 

A(lfl)«V(2>J) 

A(I^2)-V<3»JI 
A(1^3)-V(A»J) 
i-I + A 
GO TO 110 
130 CONTINIE 

C ERROR EXIT. 

WRITE (URKlNTf 2001) , kil^Z) ,1, h 

GO TO AOO 

203 WR1TE(UPRXN1»30C3) 

AGO 00 450 t-l*lUOOO 
450 REA0(UCAkD»lO01) A 

c norral exit. 

SCO REWIND UPlINCH 
C 

1001 F0RRAT(2A4>64Al ) 

20C1 FORMAT("OFRROR UNKNCWE VARIABLE "#2A1#« IN CQU'RN", I3» " ♦?"/ 
1X#2A4»64A1) 

2002 F0RNAT(1X»2A4#64A1) 

30C3 F0RnAT(**0444# cRROR. INPUT VARIABLES NOT WITHIN ALLLWED LIPITS.**# 
" SEE FR06RAP LISTING FOR LIPITS.") 

30C4 FQRHAT(**0*444« WARNING. NI IS TOO SHALL FCR STATE NOISE CASE.**) 
SToP 
END 


ccppuN lie 

CQHPUN 117 

CCPPUN lie 

CONPUN llG 
COPPUN 120 
COPPUN 121 
COHPUN 122 
COPPUN 123 
COPPUN 124 
COPPUN 125 
COPPUN 126 
COPPUN 127 
COPPUN 128 
COPPUN 129 
COPPUN 130 
COPPUN 131 
COPPUN 132 
COPPUN 133 
COPPUN 134 
COPPUN 135 
COPPUN 136 
COPPUN 137 
COPPUN 138 
COPPUN 139 
COPPUN 140 
COPPUN 141 
COPPUN 142 
COPPUN 143 
COPPUN 144 
COPPUN 145 
COPPUN 146 
COPPUN 147 
COPPUN 148 
COPPUN 149 
COPPUN 150 
COPPUN 151 
COPPUN 152 
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c 


•cjwnecK mt^tqry 
c Mnnrr tcation htstcjryi 

c nc noOTFICATIONS. 

♦ CCWDFCK AHCOH 

/AHCOM/ BC0NST(*TV) 

♦ CP»10ECK AVGCOH 

'THNiriN /AVGCOH/ ZAVGC^Z )>UAVG(»U )» EXAVC (♦EX )# ZSIG «*Z )> 

USIG(^U )#EXS1G(^EX)#ZRINH(^Z ),UM1N«(^U ) » EXMINM ( ♦EX ) , 
7niXm*7 )#UHAXH(#U ),EXHAXH( AEX) 

♦ cnMOErx «<tases 

CnMMIN /BIASES/ UQFF(^U WYOFF(#Z I 

♦ CP^O'CX BTLIN 

ClxinN /BILIN/ USEAVG»TIPVAR,Z(^Z )»U(+U )# EXTRA (♦£ X ), ONES (*6 I 
LnGT'TAL USEAVGfTIHVAP 

♦ COXOFCK '•HM 

CONiinN /COM/ NCASE#NPTT»NPT$(15)»ITMSTS(15) 

♦CPNOPCK r^MAT 

CrjWMlH /CRMAT/ AC(*X1,*X ),BC(^X1,^U )»SC(PX1,^8 )»RC(*Xlf^X )» 
CC(^ZX»PX )»^C(^Zl>^U ),HC(*Z1#PB >#ECC^Z1#^X )#FC(^XI#PX ) 
♦rOMOPCK OPTFRM 

COMM^m /OETERM/ NVAR#IMAT(*TV)fIROW(*TV)#TCOL(^TV),ILQC(^TV)» 
ACONSK^TV) 

♦CniOECK nXMMAT 

rOMMIM /DIMMAT/ a (♦Xl^^X ),B (♦XI, ♦U ),S (PXlf^B ),R (*X1,^X ), 
r (♦Zi^^x >,D (♦Z1#^U ),H (♦Z1#^B ),E (♦21#^X ) 

♦.CI3MDFCK OUMCOM 

CO»«MnN /DUMCOM/ OUM(^Xl, + X ) , 0UM2 ( ♦ X 1, *X ) , 0UM3 ( ♦Xl, ♦X ) 

♦ COMPECK OMNVEC 

COMHOM /nUMVEC/ CUMX(^X ),DUM2X(^X ) 

♦ CPMOECK PCHH 

COMMON /FCOB/ RIA(^X1,^X )#RIB(^Xi,^U ),PiS(*Xl,^B ),PI(^X1,AX ) 
^rOMOFCx PPLIST 

COMMOM /ERLIST/ BL QWUP, N XTER , E RP VE C ( 50 ) 

LOGTC*L BLOWUP 

♦ rPMPPCK CQr,M 

COMMON /FCOM/ F(^X1,^X ) 

♦ CTMOFCK *=TLES 

common /files/ UCARD,UPUNCH,UPRTNT,UCATAf LT1,UT2,UTH0UT,UWT,LPL0T 
integer UCARO,UPUNLH,UPRINT,LOATA,LTl#LT2»UTHOUTf UWTfUPLOT 

♦ COMDECX GTCOM 

COMMON /GICOM/ I TG# 01 AG6, FPE OCR, PI XG,FC1,FC2, ERRFLT# SGML S, 

GGI(PZ1,PZ ),RS0(P71,^Z ),FRS0(^Z1,^2 ),WR50(^Z ),WFRSQ(^Z ) 
LOGICAL DIAGG 
♦cpmoeck gp^ds 

COMMON /GRADS/ GR AOX ( ♦ XI , ♦NI) , GRAD Y ( ♦! 1, «N1), GR AOl I PZl, PNI) 

♦ COMOECK GPAOt 

OIMPNGIQN X0T2(^X ),X0T12(^X )#XU2(^X 1 

♦ CrMOtrK pPf^COM 

COMNON /GRDCOM/ 0K(^X1,MZ ,^KV) 

♦ CCr-OECK GRSIZE 

common /GRSIZE/ JKNMlfNK 

♦ CPMOECK MPADNG 

COMMON /HEAONG/ TI TL E ( 20 ) , ADATE, AT IM E# 

SIGLAB(2,+Z ),XLAB(2,^X ) , CON L A B( 2, ♦U ) , EX L AB ( 2, ♦EX ) 

♦ COMOECK TCONO 

COMMON /ICOND/ USERIC# VARICS,VARIC (♦X ),CX1C(^X ) 

LOGICAL USERIC,VARICS,VARIC 

♦ COMDFCK TNNAT 

COMMON /TNMAT/ ALAB,iI,JJ,IM 

♦ COMDECK TNOPT 

'“Ommon /TNOPT/ CARD,TAPE 
LO'^TCAL CARO, TAPE 

♦ COMOFCK TNORO 

CO*«MON /INORO/ NREC,ZCHAN(+Z ),UCHAN(^U ), EXC HAN ( ♦£ X ) 

TNTEGrR ZCMaN,UCHAN, EXCHAN 

♦ CPMOECK TNTFGR 

C'^mmON /INTEGR/ CT,NEAT 

♦ rOMOECK KCOM 

COMMON /KCQM/ P(^X1,^X ),KGAIN(PX1,^Z ) 

REAL KGAiN 

♦ roMDECK NAPCOH 

common /HAPrON/ PAPUK(+TV),MAPKG(^KV) 

♦C PMOECK matin 

COMMON /matin/ AV(^Xi,PX ),BV(^X1,^U ),5V(PXl,PB ),RV(^X1,PX ), 
CV(*Zi,^X ),f'V(^Zl,^U ),HV(PZi,PB ),EV(MZ1,^X ),FV(PXl,PX ), 
APRA(^X1,PX ),APRB(^X1,^U ) , A PR S ( MXl, ♦B ) , APRR (♦Xl, ♦X ), 
APRC(PZ1,PX ) ,APRD(^Z1,^U ) , A FRH( ♦Z i, ♦ B ) , APRE ( PZl, ♦X I, 
APRF(^X1,PX ),MAkD( ♦HD,7) 
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♦ COMOECt' H4TRAT 

CnHiOM /MATftAT/ )#RH(*X1,#X 

- 'T«(*Z1»#X )#OH(*Zl#*U ) 

♦ COMnfCK »»4TLAft 

C0*<inN /HATLA8/ NWaTS# LAB ( 31 ), INFL AGOl) 

REAL LAB 

♦ cr*<nErK matrix 

CnMMHM /MATRIX/ AN(AXi#AX ),BN(AXl,*C )#SN<AX1#AB )^RN(RX1,«X )p 
'‘N(AZl,*X ),0N(*Z1,RU )#HN(*Z1,AB )»EN(»Z1»*X ) 

♦CTMOFtM MAXCON 

COMinM /MAXCON/ maxhro,maxsft 

♦ <*CMnF^K MAXIM 

/maxim/ MAXX1#MAXZ1,NI,MAX1V>MAXKV 
ACrMne/*K MAXIMS 

-HMiflN /maxims/ P'AXX,MAXZ#MAXU#MAXB»LFX»LOmO 

♦ CPMDCCk MOocOM 

'‘IMMOM /MOOCQM/ UMOn 
LH';TrAL UMPf) 

♦ CTMOPCK T«^PV 

COmmHM /OBSRV/ criac(*Z1»AX )f FRIB0(RZ1#AU )#FRISH(AZ1#*B )f 
«=PICAZl,AX ) 

♦ roMOPCK imtopt 

common /OUTOPT/ FRINTX#PRlNTY»PRINTO#PlCTEM#PtNCH,TESTf PLTMAX, 

- 

LO'^T CAL PRIMTXf PRINT Yf PR INTC# PI UTcP» PUNCH# TEST# ERRTH 

♦ COMOECK RRCOM 

/FBCQM/ pmani) 

•CPMOECK PHTCHM 

COMMON /PHTCDM/ PMI(AX1#*X )»PSI(AXl#AX )# HSIB < RXl# *U )# 

- PSIS(*X1#A« ) 

♦ <*0M0FCK Rc/^rO 

<:OMMnN /RECRO/ £CFTH,T(M)#RECaPD(lCCI 
LlGTtAL EOFTH 
IMTFGCR T 
ACOMOeCM RTCC^M 

CniMHI /RICCOM/ DUMXZ(AX1,*Z ) ,DUMZX (♦ Z1»RX ), 

- 9IAP(^X1#RX I#CTG(AX1>AZ )#RIF(RX1/#X )#R1FR1F (PXl^AX ) 

♦ COMHE^K STZE 

rriMM^N /SIZE/ MX#MZ#MU#MB 

♦ COMP'^CK Sn'CQM 

CnMMPM /SOFCOM/ SnPT(*SF#7) 

♦CPMOECK SUMCOM 

COMM-IN /SUMCOM/ JKH# SUM (♦NI# AND 
♦rOMOFCK SUMSAV 

common /SUMSAV/ 0IAG0N<RNI)#APk01F(*NI )>WAPR#ITAPR 

♦ COMOPCK TAOPQS 

COMMON /TAPPOS/ ITM#REW 
LOGTCAL REW 

♦ rOMOECK TmplOTS 

OTMPHSIOH Z(*Z )#ZZtAZ )#DC(»PL)#1PLT(*PL)#VMINSC9S2)#VMAXSIRS2)# 
TIME(MTPM)# XXX(ATP*#AS2)#X(ATPA#*S )#XX(RTPR#*S ) 

EQUIVALENCE (XIUl )#XXX( 1#X) )# (XXn#l)#XXX 

c 

NTolT-ATP* 

NCN.AS 

RCOMOECX ThdaTA 

common /TOOATA/ STCn.&)#ETtU5),THIN#PRIHTl#MAXREC# 

- XBiAS(*Z I>U91AS(*U ) # EXBIAS(A£X» # ZSCALE (♦ Z )»USCALE(*U )# 
EXSCAL(REX) 

INTEGER THTN#STC»ETC 
LOGTCAL PRINT! 

♦ COMOECK XOGIRL 

COMMON /TOGIRL/ BUUNO# ERRMAX, FtJLLUNOI TER # DFaC# ITOFAC# SNOi SE 
logical FULL1#SNCISE 
♦COMOECK TOGPAO 

COMMON /T06RAD/ XT1(PX )»XT2(»X )#XH2I*X )# 

V1(#X UV2(AX >#V12(AX I#U1(*U )»L2tAU )»U12(AU )> 

- Y(AZ )»ZMY2(*Z I#ZMYFLT(*Z )#1<(AZ ) 

♦comofck toplot 

common /TDPLCT/ZPAX(»2 )#ZMIh(RZ )#XMAX(#X >,UMAX(#U ) , E XMAX ( AEX ) # 
XMIN(AX >#UMlh(*U )#EXMIN(*EX),XPL01(AX )# NUPLT# NEXPLT# 
TIMESC#RATIC#PAPER 
LOGICAL XPIOT 



• CtWO?rK VARDEFS 

•^IXt • *x 

•lAXJ • AU 
iAXQ • ♦& 

LEX • #EX 

HT - ♦NI 

• ♦TV 

mxKV • ♦‘Kv 
RAXMRO ■ ♦HO 
•lAXSFT • ♦SF 
♦roHOErK xsMHS 

Cn^HON /XSUMS/ XSUM(^X \pXZSini*% ) 

•COiO^CX FLCOSn 

Cn»<MOH /FLCONO/ CBAR , V#T HETA, PHI , A L PHA # M ACF» ♦ ARAM, C G 
RFAL »<ACH 

• COMOE^X fJPOEFS 

COMMON /GFDEFS^ GGIL AT (♦ Zl ,♦ Z )# GG IL CN (♦ 21, ♦Z ) ,FLAT(^X1#^X )# 
FLON(^XX,^X ) 

• COMDFCX GRAV 

Cni^MnH /gRAV/ 0G0T#0G0P 

• CWOECK TNFRTS 

CHNNON /inerts/ IX, I Y# IZ,IXZ,IXE,MASS, are a, CHORD, SPAN! NTCG, shift 
RFAL TX,IY,IZ, 1X2, IXE, MASS 
LOGICAL SHIFT 
^CONOECK TNSTR 

COMMON /instr/ kalf,kb,xalf,xp,xan,xax,xay,yalf,yb,yan,yax,yay, 
7ALF, 20,ZAN,ZAX,ZAY.0CGFT 
R«AL KALF,KM 

♦CON^cck lonlat 

COMMON /LONLAT/ LONG,LATR 
LOGICAL LONG,LATR 

• COMOPC*f aVCOM 

COMMON /UVCOM/ UVAR(MU ) 

^NTEGcr uvar 
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APPENDIX D 
TEST CASES 


Four test cases are provided for the MMLE3 program. These test cases were 
chosen to illustrate and check out several features of the program . This appendix 
contains descriptions of the test cases . The input cards and output listings are 
shown in supplement 2 on microfiche . Slight changes may be necessary in the 
NAMELIST format for different systems. The format shown is for CDC systems. For 
IBM systems , the dollar signs in the NAMELISTS should be changed to ampersands . 
The options and input variables illustrated in these test cases are all described in 
reference 1 , sections 3 . 3 and 4.3. 


D . 1 One-Dimensional Test Case 

The first test case uses the basic program. The system analyzed for this case 
is one-dimensional 

x(0 - Ax(t) + Bu(t) + Fn(t) x(0) = 0 

The true values of A, B, F, and G are -1, 10, 2, and 1, respectively. The 
input is a square wave with a period of 2 seconds; its value is 0 for the first second 
and 1 for the second second, repeating thereafter. Ten seconds of data at 100 samples 
per second are used. A pseudorandom number generator is used to create the state 
and measurement noise signals . Subroutine READTH is modified to compute the data 
for this test case , instead of reading a data file. This test case can be used to experi- 
ment with ways to easily implement modifications on particular computer systems . 

The case can also be run without modifying READTH by computing the data with a 
separate program . The modification instruction for READTH are shown in supple- 
ment 2 . These modifications also add subroutine GAUSSN to generate the pseudorandom 
noise. The resulting modified READTH and GAUSSN are shown in supplement 2 after 
the modification instructions . Note that no measurement noise is included in the 
first time point. Next in supplement 2 are the input cards and output listings for 
the case. 

The following options and features are used in this test case . The data channel 
numbers are specified in the NAMELIST , as the defaults are not correct for this case . 
Note that CARD, TAPE, and NREC (ref. 1, sec. 3. 3. 8 (2) and (6)) can be ignored 
since they are not used by the modified READTH (and are never used outside of 
READTH) . If (instead of being computed in READTH) the data for the test case are 
computed by a separate program and stored on a file , CARD , TAPE , and NREC will 
be relevant . ITG is set in order to turn on G determination , and the maximum total 
number of iterations is set to 10 . Note that the program converges with G fixed be- 
fore the seventh iteration , triggering the start of the G determination . Final conver- 
gence is then achieved and the program stops well before reaching the 10th iteration . 
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The TEST option is turned on for this case in order to print out dimensional matrices 
and gradients . The PRINTI option is turned on in the example so that the generated 
time history can be checked . PRINTO is turned on to provide a check on the imple- 
mentation of the time history estimation in GIRL . PRINTI and PRINTO result in large 
amounts of output (particularly PRINTI); thus, this case should probably be run 
first without them. The case can be rerun with them turned on if necessary for de- 
bugging. The RELAB option is used to read in more meaningful labels than the 
defaults . The time history plot option is on by default . A plot is requested of the 
corrected state estimate . NEXPLT is used to request plots of the true state , state 
noise , and measurement noise , which are carried as the first three extra signals for 
this case. Setting NUPLT to 1 does not affect the resulting plot for this case, but 
saves the computer time otherwise required to notice that the second to fourth controls 
are identically zero and thus need not be plotted . 

Starting values are read in for the AN , CN , F , and GGI matrices . Both full and 
diagonal input formats are used (there is little difference for 1X1 matrices) . A 
starting BN is not read in , so it is 0 , along with all of the other unmentioned matrices 
except for RN . The AV , BV , and FV matrices are read in to specify the unknowns . 
GGI is treated separately and is specified to be unknown by ITG in the NAMELIST. 

The values of MX and MZ are set to 1 from the sizes read in for the AN and GGI mat- 
rices . The BN matrix was not read in, but MU is set to 1 because of the unknown 
derivative of the first control specified by BV . The value of MB is set to NCASE 
which is 1; even though S and H are both 0 for this case, MB must be defined, and 0 
is not an allowable value . 

The final estimates , Cramer-Rao bounds , and true values are found in the 
following table . 



Estimate 

Cramer-Rao 

bound 

True value 

AN 

-0.8897 

0.20 

-1.0 

BN 

9.619 

1.5 

10.0 

FN 

2.071 

0.21 

2.0 


The measurement noise covariance matrix, is not directly estimated, but 

can be computed from the residual power estimate , GGI ^ , and the estimate of the 
prediction error power , P (the Riccati covariance matrix) . 

^ = Vg^ - P 
= Vggi"^ - P 
= >.8176"^ - .2184 
= 1.0023 
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This compares to the true value of 1 . A Cramer-Rao bound is not computed 
for 


The time history plots from this case are shown below . 



105 


D.2 


D.2 Longitudinal Test Cases 

The second and third test cases are longitudinal cases using the standard 
aircraft routines. Both use actual flight data from a T-37 aircraft (ref. 10) . These 
two cases are set up as a single job , but can be run separately . 

A common type of update to the standard aircraft routines is illustrated in these 
check cases: a modification to automatically compute weights and inertias . Some- 

times fuel weights or other quantities recorded as extra signals are used for such 
computations . In the test case here , a table of fuel weights and times obtained from 
pilot lap notes is read in subroutine ONCE. Also read in are tables to obtain total 
weight and inertias as a function of fuel weight . Subroutine AVERAG then uses 
linear interpolation on these tables to automatically obtain the total weight and inertias 
for each maneuver . The modification instructions for this update , followed by the 
resulting modified program listings , are shown in supplement 2 . The test cases can 
be run without modifying the program by entering the weights and inertias from the 
output listing into NAMELIST USER. 

Next on supplement 2 are the input cards and output listings for the test cases. 

A simple predicted- derivative data set is used which has a constant value for each 
derivative. One lateral- directional derivative, , is included to illustrate that it 
is ignored for these longitudinal cases. P 

Test case 2 is an unusual maneuver, designed for estimating C (ref. 11) . 

a 

This maneuver requires several of the program's more sophisticated features. 

Because of the significant variations of q, V, 0, and cp during the maneuver, the 

time-varying option must be used. Since the aircraft does a complete 360° roll 

during this longitudinal maneuver, the lateral-directional cross-coupling terms are 

quite important , particularly in the a and 0 equations . MZ is set to 4 to eliminate the 

a observation equation , using the simplified low a longitudinal equations (compare 
oc 

the lists of unknowns in the output of the second and third test cases) . In order to 
reasonably match a for this maneuver , thrust and drag must be treated separately 

because of the q variation. (Such treatment can be made, but is not included in the 

test case.) The default automatic scaling for altitude (extra signal 7) is overridden 

to obtain a more sensitive scale , since the default would include 0 in the scale . 

FREQCR is set in order to obtain filtered residual powers and use them to adjust the 

Cramer-Rao bounds . An RV matrix is read in to define C (RN„ ,) to be unknown. 

m« 2,1 
a 

The time history follows on cards . The record length of the input time history is set 
to 23 in order to use one less card per time point than with the default, because the 
last two channels are not used for this case . 

Convergence is rapid and monotone . The first iteration changes the linear 

unknowns only and lowers the cost functional by 2^ orders of magnitude . The second 

iteration changes all of the unknowns and reduces the cost functional by a factor of 4 . 
The solution has then been essentially reached; the remaining iterations just add more 
significant digits. Note that the filtered error sum and log determinant are slightly 
better (in the fourth and fifth places) after iteration 3 than at the final value . This 
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is not unusual , because these are not the quantities minimized (the unfiltered 
weighted error sum is minimized) . Needless to say , these quantities should be 
expected to be at least near their minima as they are on this case . The fit is 
shown on the following pages . 
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Test case 3 is a set of elevator pulses typical of the data normally obtained for 
longitudinal stability and control derivative estimation . The time- varying option is 
not needed for this case. The more complicated longitudinal equations including a 

tX, 

are used in the sample run (although the case will run quite well and use less com- 
puter time if MZ is set to 4 to eliminate a ) . The convergence is excellent and similar 

to the previous case. Several observations can be made about the fit shown below. 



-10 •• 

0 1 2 3 4 6 6 7 

TIME 

MMLE3 TEST CASE 3. T37 FLT 190 CASE 20-21. UP AND DOHM ELEVATOR PULSES- 


11 . 34 . 9 . 23 . 



D.3 


The resolution on 0 is relatively poor for this maneuver , but does not appear to 
have caused any problems (the rounded corners on the bit jumps are due to digital 
filtering done after the flight) . The a match shows some significant discrepancies . 

The discrepancies are strongly correlated with a, and not with q or 6^. This sug- 
gests a significant nonlinearity in the versus a curve. Since a ranges from -1° 

to 5°, it would not be surprising to find the linear derivative inadequate. A 2 

a 

term, either replacing or in addition to the term, might give a better model. 

« 2 

Nonlinear terms such as C„ „ can be implemented in MMLE3 by forming a as an extra 

control . This task is appropriate for subroutine THMOD or can be done by a separate 
program to add the signal to the data file (of course , the a used for this purpose 
should be corrected to the center of gravity and for any upwash) . Running this case 
with a 2 term is left as a relatively simple exercise in using MMLE3. 


D . 3 Lateral-Directional Test Case 

The fourth test case is a lateral-directional case using the standard aircraft 
routines . This case consists of actual flight data from an oblique wing aircraft 
(ref. 12); the wing is not skewed during this maneuver. The data are typical of 
those obtained for lateral-directional stability and control derivative estimation . 

There are no program modifications for this test case. The input cards and 
output listing are in supplement 2 following the previous test cases . This case is 
run without a predicted-derivative data set to illustrate that option. Vehicle geometry 
and instrument positions are read in NAMELIST USER. The weight and inertias are 
also read in the NAMELIST for this case , in contrast to the previous two test cases , 
which contain an update to compute them . This case is run using metric units . The 
engine revolutions per minute (extra signal 11) is specified to be found on channel 10 
of the input data instead of the default channel 26 . Since no predicted-derivative 
data set was used , a starting AN matrix must be read in to provide reasonable starting 
estimates . Note that the AN read in is not square. The program accepts this input, 
but later changes the dimensions used to be consistent. No BN matrix is read in, so 
the starting estimates of the control derivatives are all 0 . BV is read in order to 
override the default that includes as an unknown. The BV matrix read in also 

a 

includes values in the fifth row to illustrate that they are ignored. The number of 
rows of BN (and thus BV) is forced internally to equal MX (4 for this case) and any 
entries outside of this range are ignored . 

Convergence is rapid and uneventful . The resulting fit is shown below . 


Ill 
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